0.0.1 Preliminaries about Spatial Data

#Spatial Data 在R上的图形展示
#https://zhuanlan.zhihu.com/p/27547697
#R code for Moran's I and Geary's C
#clean up
rm(list=ls())

#preferences
options(scipen=8)

#load required libraries
library(spdep) 
## Loading required package: sp
## Loading required package: Matrix
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(maps)
library(maptools)
## Checking rgeos availability: TRUE
library(foreign)
library(rgdal)
## rgdal: version: 1.3-6, (SVN revision 773)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/gdal
##  GDAL binary built with GEOS: FALSE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/proj
##  Linking to sp version: 1.3-1
#Load Data
imm.data.0<-read.dta("stateImmig0511.dta")
#fix(imm.data.0)
#names(imm.data.0)

#subset
imm.data<-subset(imm.data.0,State!="ALASKA" & State!="HAWAII")

#Estimate a Linear Model
mod.immigrant<-lm(immig0511~pubIdeolCCES, data=imm.data); summary(mod.immigrant)
## 
## Call:
## lm(formula = immig0511 ~ pubIdeolCCES, data = imm.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.33009 -0.55823  0.04261  0.48461  1.38787 
## 
## Coefficients:
##              Estimate Std. Error t value     Pr(>|t|)    
## (Intercept)   0.93593    0.16676   5.613 0.0000010985 ***
## pubIdeolCCES  0.06738    0.01058   6.370 0.0000000806 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6671 on 46 degrees of freedom
## Multiple R-squared:  0.4687, Adjusted R-squared:  0.4571 
## F-statistic: 40.58 on 1 and 46 DF,  p-value: 0.00000008058
#call data from "maps", create neighbor matrix
us.states<-map("state", fill=TRUE, plot=TRUE, resolution=0) #too cluttered, projection=???

IDs <- sapply(strsplit(us.states$names, ":"), function(x) x[1])
us.map.0<-map2SpatialPolygons(us.states, IDs=IDs)

#eliminate Washington, DC
us.map<-us.map.0[-8]

#create neighbor list from state map
state.nb<-poly2nb(us.map)
state.coords<-coordinates(us.map)
plot(state.nb, coords=state.coords)

#state by state list, see how it compares to Melanie Wall's list: http://www.biostat.umn.edu/~brad/data/contig-lower48.dat
for(i in 1:48){print(state.nb[i])}
## [[1]]
## [1]  8  9 22 40
## 
## [[1]]
## [1]  4  5 26 29 42
## 
## [[1]]
## [1] 16 22 23 34 40 41
## 
## [[1]]
## [1]  2 26 35
## 
## [[1]]
## [1]  2 14 25 29 34 42 48
## 
## [[1]]
## [1] 19 30 37
## 
## [[1]]
## [1] 18 28 36
## 
## [[1]]
## [1] 1 9
## 
## [[1]]
## [1]  1  8 31 38 40
## 
## [[1]]
## [1] 24 26 35 42 45 48
## 
## [[1]]
## [1] 12 13 15 23 47
## 
## [[1]]
## [1] 11 15 20 33
## 
## [[1]]
## [1] 11 21 23 25 39 47
## 
## [[1]]
## [1]  5 23 25 34
## 
## [[1]]
## [1] 11 12 23 33 40 44 46
## 
## [[1]]
## [1]  3 22 41
## 
## [[1]]
## [1] 27
## 
## [[1]]
## [1]  7 36 44 46
## 
## [[1]]
## [1]  6 27 30 37 43
## 
## [[1]]
## [1] 12 33 47
## 
## [[1]]
## [1] 13 32 39 47
## 
## [[1]]
## [1]  1  3 16 40
## 
## [[1]]
## [1]  3 11 13 14 15 25 34 40
## 
## [[1]]
## [1] 10 32 39 48
## 
## [[1]]
## [1]  5 13 14 23 39 48
## 
## [[1]]
## [1]  2  4 10 35 42
## 
## [[1]]
## [1] 17 19 43
## 
## [[1]]
## [1]  7 30 36
## 
## [[1]]
## [1]  2  5 34 41 42
## 
## [[1]]
## [1]  6 19 28 36 43
## 
## [[1]]
## [1]  9 38 40 44
## 
## [[1]]
## [1] 21 24 39
## 
## [[1]]
## [1] 12 15 20 36 46
## 
## [[1]]
## [1]  3  5 14 23 29 41
## 
## [[1]]
## [1]  4 10 26 45
## 
## [[1]]
## [1]  7 18 28 30 33 46
## 
## [[1]]
## [1]  6 19
## 
## [[1]]
## [1]  9 31
## 
## [[1]]
## [1] 13 21 24 25 32 48
## 
## [[1]]
## [1]  1  3  9 15 22 23 31 44
## 
## [[1]]
## [1]  3 16 29 34
## 
## [[1]]
## [1]  2  5 10 26 29 48
## 
## [[1]]
## [1] 19 27 30
## 
## [[1]]
## [1] 15 18 31 40 46
## 
## [[1]]
## [1] 10 35
## 
## [[1]]
## [1] 15 18 33 36 44
## 
## [[1]]
## [1] 11 13 20 21
## 
## [[1]]
## [1]  5 10 24 25 39 42
### MORAN'S I ###
#First, is the dependent variable spatially correlated on its own?
moran.test(imm.data$immig0511,nb2listw(state.nb,style='W'))
## 
##  Moran I test under randomisation
## 
## data:  imm.data$immig0511  
## weights: nb2listw(state.nb, style = "W")    
## 
## Moran I statistic standard deviate = 3.7071, p-value = 0.0001048
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.34159785       -0.02127660        0.00958191
#Moran's I on residuals
moran.test(resid(mod.immigrant),nb2listw(state.nb,style='W'))
## 
##  Moran I test under randomisation
## 
## data:  resid(mod.immigrant)  
## weights: nb2listw(state.nb, style = "W")    
## 
## Moran I statistic standard deviate = 0.42017, p-value = 0.3372
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.019850914      -0.021276596       0.009581157
# Moran's I test for residuals from a regression model
#   based on weighted residuals
lm.morantest(mod.immigrant,nb2listw(state.nb,style='W'))
## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = immig0511 ~ pubIdeolCCES, data = imm.data)
## weights: nb2listw(state.nb, style = "W")
## 
## Moran I statistic standard deviate = 0.54009, p-value = 0.2946
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.019850914     -0.031418848      0.009011362
### GEARY'S C ###
#Raw dependent variable
geary.test(imm.data$immig0511,nb2listw(state.nb,style='W'))
## 
##  Geary C test under randomisation
## 
## data:  imm.data$immig0511 
## weights: nb2listw(state.nb, style = "W") 
## 
## Geary C statistic standard deviate = 4.0725, p-value = 0.00002326
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.59254949        1.00000000        0.01001005
#Geary's C on residuals
geary.test(resid(mod.immigrant),nb2listw(state.nb,style='W'))
## 
##  Geary C test under randomisation
## 
## data:  resid(mod.immigrant) 
## weights: nb2listw(state.nb, style = "W") 
## 
## Geary C statistic standard deviate = 0.77845, p-value = 0.2182
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.92211018        1.00000000        0.01001147
#For homework:
#A tip for setting-up your neighbor matrix for a state using a map from the "maps" library
iowa<-map("county", "iowa", fill=TRUE, plot=TRUE, resolution=0) 

IDs <- sapply(strsplit(iowa$names, ","), function(x) x[2])
ia.map.0<-map2SpatialPolygons(iowa, IDs=IDs) 

#Also for homework:
#Type ?cos for more information on trigonometric functions

0.0.2 Regression with Spatially Lagged Dependent Variables

###WORKING EXAMPLE FROM CHAPTER 2 OF WARD & GLEDITSCH 2008###
###FUTURE USE OF THIS CODE REQUIRES CITATION OF THEIR ORIGINAL WORK###

#Load Data
# sldv is the data.frame
# mdd2 is the minimum distance data.frame
source("https://spia.uga.edu/faculty_pages/monogan/teaching/spatial/chapter1data.R")
#fix(sldv)
head(sldv)
##   polygon.number  place.name tla country.name population kilometer2
## 1              3 Afghanistan AFG  Afghanistan   17250390  641869.19
## 2              6     Albania ALB      Albania    3416945   28754.50
## 3              4     Algeria ALG      Algeria   27459230 2320972.00
## 4              9      Angola ANG       Angola   11527260 1252421.00
## 5             11   Argentina ARG    Argentina   33796870 2781013.00
## 6              7     Armenia ARM      Armenia    3377228   29872.46
##   idnumber     gdp.2002 democracy
## 1        3  20000000000       -10
## 2        6   4840000000         5
## 3        4  55900000000        -3
## 4        9  11200000000        -3
## 5       11 102000000000         8
## 6        7   2370000000         5
#fix(mdd2)
head(mdd2)
##        ida idb year mindist
## 114585 AFG CHN 2002       0
## 114586 AFG IND 2002     241
## 114587 AFG IRN 2002       0
## 114588 AFG KYR 2002     116
## 114589 AFG KZK 2002     331
## 114590 AFG OMA 2002     730
# Create a neighbor minimum distance matrix
# from the dyadic minimum distance data
# (Note that the cshapes library provides a 
# convienent way to create a matrix directly)
mddmat <- matrix(9999,ncol=dim(sldv)[1],nrow=dim(sldv)[1]) #We choose 9999km as the default because it's a big number that indicates "not a neighbor."
dimnames(mddmat) <- list(c(sldv$tla),c(sldv$tla))
for(i in 1:dim(mdd2)[1]){
    mddmat[mdd2$ida[i],mdd2$idb[i]] <-  mdd2$mindist[i] #super cool stuff: call text entries to set matrix cells
}
mddmat[1:12,1:12]
##      AFG  ALB  ALG  ANG  ARG  ARM  AUL  AUS  AZE  BAH  BEL  BEN
## AFG 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999
## ALB 9999 9999 9999 9999 9999 9999 9999  554 9999 9999 9999 9999
## ALG 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999  763
## ANG 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999
## ARG 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999
## ARM 9999 9999 9999 9999 9999 9999 9999 9999    0 9999 9999 9999
## AUL 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999
## AUS 9999  554 9999 9999 9999 9999 9999 9999 9999 9999  369 9999
## AZE 9999 9999 9999 9999 9999    0 9999 9999 9999 9999 9999 9999
## BAH 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999
## BEL 9999 9999 9999 9999 9999 9999 9999  369 9999 9999 9999 9999
## BEN 9999 9999  763 9999 9999 9999 9999 9999 9999 9999 9999 9999
# Create a binary matrix with a 200 km threshold
m200mat <- mddmat
m200mat[m200mat<=200] <- 1
m200mat[m200mat>200] <- 0
m200mat[1:12,1:12]
##     AFG ALB ALG ANG ARG ARM AUL AUS AZE BAH BEL BEN
## AFG   0   0   0   0   0   0   0   0   0   0   0   0
## ALB   0   0   0   0   0   0   0   0   0   0   0   0
## ALG   0   0   0   0   0   0   0   0   0   0   0   0
## ANG   0   0   0   0   0   0   0   0   0   0   0   0
## ARG   0   0   0   0   0   0   0   0   0   0   0   0
## ARM   0   0   0   0   0   0   0   0   1   0   0   0
## AUL   0   0   0   0   0   0   0   0   0   0   0   0
## AUS   0   0   0   0   0   0   0   0   0   0   0   0
## AZE   0   0   0   0   0   1   0   0   0   0   0   0
## BAH   0   0   0   0   0   0   0   0   0   0   0   0
## BEL   0   0   0   0   0   0   0   0   0   0   0   0
## BEN   0   0   0   0   0   0   0   0   0   0   0   0
# Create a row standardized weights matrix
wmat <- m200mat/apply(m200mat,1,sum)
wmat[1:12,1:12]
##     AFG ALB ALG ANG ARG       ARM AUL AUS AZE BAH BEL BEN
## AFG   0   0   0   0   0 0.0000000   0   0 0.0   0   0   0
## ALB   0   0   0   0   0 0.0000000   0   0 0.0   0   0   0
## ALG   0   0   0   0   0 0.0000000   0   0 0.0   0   0   0
## ANG   0   0   0   0   0 0.0000000   0   0 0.0   0   0   0
## ARG   0   0   0   0   0 0.0000000   0   0 0.0   0   0   0
## ARM   0   0   0   0   0 0.0000000   0   0 0.2   0   0   0
## AUL   0   0   0   0   0 0.0000000   0   0 0.0   0   0   0
## AUS   0   0   0   0   0 0.0000000   0   0 0.0   0   0   0
## AZE   0   0   0   0   0 0.1428571   0   0 0.0   0   0   0
## BAH   0   0   0   0   0 0.0000000   0   0 0.0   0   0   0
## BEL   0   0   0   0   0 0.0000000   0   0 0.0   0   0   0
## BEN   0   0   0   0   0 0.0000000   0   0 0.0   0   0   0
# calculate the spatial lag of democracy
democracy.spatial.lag <- as.vector(wmat%*%sldv$democracy)
democracy.spatial.lag[1:12]
##  [1]  -4.5714286   7.0000000  -0.4285714   2.2000000   8.6000000
##  [6]   3.0000000   8.5000000   8.5454545   1.7142857 -10.0000000
## [11]   9.8000000   2.4000000
# Create a weights list object (a list object is required for spdep)
cmat.lw <- spdep::mat2listw(m200mat, row.names=row.names(m200mat))

# Extract a neighbor list
cmat.nb <- cmat.lw$neighbours




### REALLY NOT RECOMMENDED: OLS WITH NO SPATIAL INFORMATION ###
stupid.model<-lm(democracy~log(gdp.2002/population), data=sldv); summary(stupid.model)
## 
## Call:
## lm(formula = democracy ~ log(gdp.2002/population), data = sldv)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.945  -4.942   2.557   4.544   9.251 
## 
## Coefficients:
##                          Estimate Std. Error t value    Pr(>|t|)    
## (Intercept)               -9.6874     2.4258  -3.994 0.000099996 ***
## log(gdp.2002/population)   1.6780     0.3128   5.364 0.000000289 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.285 on 156 degrees of freedom
## Multiple R-squared:  0.1557, Adjusted R-squared:  0.1503 
## F-statistic: 28.77 on 1 and 156 DF,  p-value: 0.0000002892
moran.test(resid(stupid.model),nb2listw(cmat.nb))
## 
##  Moran I test under randomisation
## 
## data:  resid(stupid.model)  
## weights: nb2listw(cmat.nb)    
## 
## Moran I statistic standard deviate = 7.7664, p-value = 4.038e-15
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.404450410      -0.006369427       0.002798125
### NOT RECOMMENDED: OLS WITH A SPATIAL LAG ###
sldv$dem.spat.lag<-democracy.spatial.lag
silly.model<-lm(democracy~log(gdp.2002/population)+dem.spat.lag, 
                data=sldv); summary(silly.model)
## 
## Call:
## lm(formula = democracy ~ log(gdp.2002/population) + dem.spat.lag, 
##     data = sldv)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.2447  -3.5407   0.6011   2.8480  11.7758 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              -4.97495    2.07116  -2.402  0.01749 *  
## log(gdp.2002/population)  0.75927    0.27872   2.724  0.00719 ** 
## dem.spat.lag              0.76171    0.08802   8.653    6e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.177 on 155 degrees of freedom
## Multiple R-squared:  0.4307, Adjusted R-squared:  0.4234 
## F-statistic: 58.64 on 2 and 155 DF,  p-value: < 2.2e-16
moran.test(resid(silly.model),nb2listw(cmat.nb))
## 
##  Moran I test under randomisation
## 
## data:  resid(silly.model)  
## weights: nb2listw(cmat.nb)    
## 
## Moran I statistic standard deviate = -3.1255, p-value = 0.9991
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.171459851      -0.006369427       0.002789921
#But do not chase all that glitters... How are we doing Gauss-Markov wise? Poorly.

### BETTER CHOICE: MLE WITH A SPATIAL LAG ###
sldv.fit <- spdep::lagsarlm(democracy ~ log(gdp.2002/population), 
                     data=sldv, nb2listw(cmat.nb), 
                     method="eigen", quiet=FALSE); summary(sldv.fit)
## 
## Spatial lag model
## Jacobian calculated using neighbourhood matrix eigenvalues
## Computing eigenvalues ...
## 
## rho:  -0.4815483     function value:  -551.2436 
## rho:  0.0843528  function value:  -508.153 
## rho:  0.4340989  function value:  -492.6041 
## rho:  0.6502539  function value:  -491.8775 
## rho:  0.565321   function value:  -491.0999 
## rho:  0.559797   function value:  -491.1005 
## rho:  0.5631131  function value:  -491.0994 
## rho:  0.5631432  function value:  -491.0994 
## rho:  0.5631453  function value:  -491.0994 
## rho:  0.5631453  function value:  -491.0994 
## rho:  0.5631453  function value:  -491.0994 
## rho:  0.5631453  function value:  -491.0994
## 
## Call:spdep::lagsarlm(formula = democracy ~ log(gdp.2002/population), 
##     data = sldv, listw = nb2listw(cmat.nb), method = "eigen", 
##     quiet = FALSE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.2640  -4.1620   1.2861   3.5655  11.1175 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                          Estimate Std. Error z value  Pr(>|z|)
## (Intercept)              -6.20342    2.08233 -2.9791 0.0028912
## log(gdp.2002/population)  0.99877    0.27826  3.5894 0.0003315
## 
## Rho: 0.56315, LR test value: 45.035, p-value: 0.000000000019358
## Asymptotic standard error: 0.07576
##     z-value: 7.4333, p-value: 1.0592e-13
## Wald statistic: 55.254, p-value: 1.0592e-13
## 
## Log likelihood: -491.0994 for lag model
## ML residual variance (sigma squared): 27.161, (sigma: 5.2116)
## Number of observations: 158 
## Number of parameters estimated: 4 
## AIC: 990.2, (AIC for lm: 1033.2)
## LM test for residual autocorrelation
## test value: 2.1028, p-value: 0.14703
spdep::moran.test(resid(sldv.fit),nb2listw(cmat.nb))
## 
##  Moran I test under randomisation
## 
## data:  resid(sldv.fit)  
## weights: nb2listw(cmat.nb)    
## 
## Moran I statistic standard deviate = -0.45541, p-value = 0.6756
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.030459980      -0.006369427       0.002798233
### CALCULATING EQUILIBRIUM EFFECTS ###
# Code to calculate equilibrium effect of changes in GDP per capita
# Create vector to store the estimate for each states
ee.est <- rep(NA,dim(sldv)[1])

# Assign the country name labels
names(ee.est) <- sldv$tla

# Create a null vector to use in loop
svec <- rep(0,dim(sldv)[1])

# Create a N x N identity matrix
eye <- matrix(0,nrow=dim(sldv)[1],ncol=dim(sldv)[1])
diag(eye) <- 1

# Loop over 1:n states and store effect of change in
# each state i in ee.est[i]
for(i in 1:length(ee.est)){
    cvec <- svec
    cvec[i] <- 1
    res <- solve(eye - 0.56315 * wmat) %*% cvec * 0.99877
    ee.est[i] <- res[i]
}

## Results by country: How much would democracy increase in a country if log(gdp.2002/population) increased by 1 in that country? ##
ee.est
##      AFG      ALB      ALG      ANG      ARG      ARM      AUL      AUS 
## 1.067492 1.064252 1.104560 1.071895 1.141127 1.046587 1.121907 1.068377 
##      AZE      BAH      BEL      BEN      BFO      BHM      BHU      BLR 
## 1.054909 1.063698 1.089260 1.085159 1.093464 1.103532 1.077306 1.065284 
##      BLZ      BNG      BOL      BOS      BOT      BRA      BUI      BUL 
## 1.120205 1.084170 1.121271 1.057845 1.086653 1.160317 1.077169 1.060156 
##      CAM      CAN      CAO      CDI      CEN      CHA      CHL      CHN 
## 1.076181 1.071441 1.097271 1.089857 1.059433 1.061780 1.090393 1.104231 
##      COL      CON      COS      CRO      CUB      CYP      CZR      DEN 
## 1.122435 1.080874 1.114834 1.057845 1.151169 1.055224 1.043784 1.059586 
##      DJI      DOM      DRC      DRV      ECU      EGY      EQG      ERI 
## 1.094217 1.095957 1.086186 1.100169 1.074091 1.059802 1.081722 1.092907 
##      EST      ETH      FIN      FJI      FRN      GAB      GAM      GFR 
## 1.076870 1.097106 1.067623 1.130221 1.093156 1.089852 1.107729 1.087953 
##      GHA      GNB      GRC      GRG      GUA      GUI      GUY      HAI 
## 1.084570 1.107729 1.056691 1.049064 1.120205 1.154450 1.122284 1.225806 
##      HON      HUN      IND      INS      IRE      IRN      IRQ      ISR 
## 1.126266 1.063546 1.131937 1.151724 1.090164 1.080878 1.064917 1.073833 
##      ITA      JAM      JOR      JPN      KEN      KUW      KYR      KZK 
## 1.070141 1.101928 1.079516 1.108985 1.073287 1.043256 1.059848 1.070691 
##      LAO      LAT      LBR      LEB      LES      LIB      LIT      LUX 
## 1.084320 1.069949 1.120316 1.077060 1.060203 1.070449 1.064518 1.069460 
##      MAA      MAC      MAL      MAW      MEX      MLD      MLI      MON 
## 1.087258 1.060525 1.094080 1.066915 1.108328 1.034501 1.080262 1.028817 
##      MOR      MYA      MZM      NAM      NEP      NIC      NIG      NIR 
## 1.090795 1.071990 1.105845 1.086653 1.077306 1.139634 1.086168 1.078120 
##      NOR      NTH      OMA      PAK      PAN      PAR      PER      PHI 
## 1.074484 1.089260 1.064920 1.046347 1.129586 1.080619 1.119956 1.104765 
##      PNG      POL      POR      PRK      QAT      ROK      RUM      RUS 
## 1.243478 1.065041 1.090252 1.073413 1.092376 1.121669 1.068240 1.089609 
##      RWA      SAF      SAL      SAU      SEN      SIE      SLO      SLV 
## 1.067805 1.162269 1.114806 1.116874 1.118326 1.099129 1.055019 1.052143 
##      SOM      SPN      SRI      SUD      SUR      SWA      SWD      SWZ 
## 1.098215 1.153850 1.050053 1.074267 1.090706 1.067901 1.080689 1.039917 
##      SYR      TAJ      TAW      TAZ      THI      TKM      TOG      TRI 
## 1.081751 1.073590 1.087880 1.088571 1.094611 1.056263 1.085159 1.090907 
##      TUN      TUR      UAE      UGA      UKG      UKR      URU      USA 
## 1.050776 1.071617 1.076321 1.077890 1.152739 1.065221 1.064233 1.145732 
##      UZB      VEN      YEM      YUG      ZAM      ZIM 
## 1.066881 1.162103 1.101142 1.069267 1.092730 1.096897
hist(ee.est)

## If log(gdp.2002/population) increased by 1 in Russia, what would happen to democracy in other states (observation 120)? ##
cvec <- rep(0,dim(sldv)[1])
cvec[120] <- 1

# Store estimates for impact of change in Russia in rus.est
eye <- matrix(0,nrow=dim(sldv)[1],ncol=dim(sldv)[1])
diag(eye) <- 1
rus.est <- solve(eye - 0.56315 * wmat) %*% cvec*0.99877

# Find ten highest values of rus.est vector
rus.est <- round(rus.est,3)
rus.est <- data.frame(sldv$tla,rus.est)
rus.est[rev(order(rus.est$rus.est)),][1:10,]
##     sldv.tla rus.est
## RUS      RUS   1.090
## PRK      PRK   0.242
## JPN      JPN   0.241
## MON      MON   0.240
## FIN      FIN   0.222
## EST      EST   0.212
## NOR      NOR   0.202
## LIT      LIT   0.198
## LAT      LAT   0.197
## ARM      ARM   0.177
## What if China had a democracy score of 10? How would other states change in democracy? ##
# China is observation 32
cvec <- rep(0,dim(sldv)[1])
cvec[32] <- 10

# Store estimates of change in China in chn.est
chn.est <- c(cbind (0, 0, wmat%*%cvec) %*% c(summary(sldv.fit)$Coef[,1],summary(sldv.fit)$rho))
chn.est <- round(chn.est,3)

# Find all states where non-zero impact
chn.est <- data.frame(sldv$tla,chn.est)
chn.est.2<-subset(chn.est,chn.est>0)
chn.est.2[rev(order(chn.est.2$chn.est)),]
##     sldv.tla chn.est
## 139      TAW   1.877
## 116      PRK   1.877
## 96       MON   1.877
## 101      NEP   1.408
## 15       BHU   1.408
## 108      PAK   1.126
## 81       LAO   1.126
## 79       KYR   1.126
## 18       BNG   1.126
## 153      UZB   0.939
## 141      THI   0.939
## 98       MYA   0.939
## 138      TAJ   0.804
## 67       IND   0.804
## 44       DRV   0.804
## 1        AFG   0.804
## 80       KZK   0.704
## 120      RUS   0.282
#### TIPS FOR HOMEWORK ####
# (1) #
#You have to load "chapter1data.R" again to get the "mdd2" information. We have not use for "sldv", though. After running "chapter1data.R", to avoid confusion, I'd type:
rm(sldv)

# (2) #
neighborhood.0<-read.csv("neighbor2002subset.csv")
#After you load the "neighbor2002subset.csv" dataset (which I named "neighborhood.0"), run the following line of code:
modeled<-!is.na(neighborhood.0$polity2) & 
  !is.na(neighborhood.0$lngdppercap) & 
  !is.na(neighborhood.0$physint)
#Create a second copy of your data subsetting on the conditions of "modeled".

# (3) #
#When defnining "mddmat" use code like this:
included<-c(as.character(neighborhood.0$tla))
mddmat <- matrix(9999,ncol=dim(neighborhood.0)[1],
                 nrow=dim(neighborhood.0)[1]) 
dimnames(mddmat) <- list(included,included)
for(i in 1:dim(mdd2)[1]){
    if(is.element(mdd2$ida[i],included) & is.element(mdd2$idb[i],included)){
        mddmat[mdd2$ida[i],mdd2$idb[i]] <-  mdd2$mindist[i] #super cool stuff: call text entries to set matrix cells
        }
}
#Here, "neighborhood" is my subsetted data frame. Convoluted code, I know, but I'll try to explain.

0.0.3 The Spatial Error Model

#clean up
rm(list=ls())

#preferences
options(scipen=8)

#load required libraries
library(spdep) 
library(maps)
library(foreign)
library(rgdal)


#Load Data
imm.data.0<-read.dta("stateImmig0511.dta")
#names(imm.data.0)

#subset
imm.data<-subset(imm.data.0,State!="ALASKA" & State!="HAWAII")

#call data from "maps", create neighbor matrix
us.states<-map("state", fill=TRUE, plot=TRUE, resolution=0) 

IDs <- sapply(strsplit(us.states$names, ":"), function(x) x[1])
us.map.0<-map2SpatialPolygons(us.states, IDs=IDs)

#eliminate Washington, DC
us.map<-us.map.0[-8]

#create neighbor list from state map
state.nb<-poly2nb(us.map)
state.coords<-coordinates(us.map)
plot(state.nb, coords=state.coords)

###estimate a simple linear model###
mod.immigrant<-lm(immig0511~pubIdeolCCES, data=imm.data); summary(mod.immigrant)
## 
## Call:
## lm(formula = immig0511 ~ pubIdeolCCES, data = imm.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.33009 -0.55823  0.04261  0.48461  1.38787 
## 
## Coefficients:
##              Estimate Std. Error t value     Pr(>|t|)    
## (Intercept)   0.93593    0.16676   5.613 0.0000010985 ***
## pubIdeolCCES  0.06738    0.01058   6.370 0.0000000806 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6671 on 46 degrees of freedom
## Multiple R-squared:  0.4687, Adjusted R-squared:  0.4571 
## F-statistic: 40.58 on 1 and 46 DF,  p-value: 0.00000008058
###Geary's C on residuals###
geary.test(resid(mod.immigrant),nb2listw(state.nb,style='C'))
## 
##  Geary C test under randomisation
## 
## data:  resid(mod.immigrant) 
## weights: nb2listw(state.nb, style = "C") 
## 
## Geary C statistic standard deviate = 1.5279, p-value = 0.06327
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##         0.8346617         1.0000000         0.0117107
###plot residuals###
X <- c(mod.immigrant$residuals[1:20],mod.immigrant$residuals[20],mod.immigrant$residuals[21:48])
n.col <- 10
br <- c(-999,quantile(X,c(.1,.2,.3,.4,.5,.6,.7,.8,.9)), 999)
#shading<-gray((n.col-1):0/(n.col-1))
shading<- colorRampPalette(c("red","white","blue"))(10)
#shading<- colorRampPalette(c("red","blue"))(10)

X.grp<-findInterval(X, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.shad<-shading[X.grp]
states.to.map<-c("alabama","arizona","arkansas","california",
                 "colorado","connecticut","delaware","florida",
                 "georgia","idaho","illinois","indiana","iowa",
                 "kansas","kentucky","louisiana","maine",
                 "maryland","massachusetts:main","michigan:north",
                 "michigan:south","minnesota","mississippi",
                 "missouri","montana","nebraska","nevada",
                 "new hampshire","new jersey","new mexico",
                 "new york:main","north carolina:main","north dakota",
                 "ohio","oklahoma","oregon","pennsylvania",
                 "rhode island","south carolina","south dakota","tennessee",
                 "texas","utah","vermont","virginia:main",
                 "washington:main","west virginia","wisconsin","wyoming")
map(database="state", region=states.to.map, fill=T, 
    plot=T, interior=T, col=X.shad, exact=T)

#Fit spatial error model
spat.err.immig<-errorsarlm(immig0511~pubIdeolCCES, data=imm.data,
                           nb2listw(state.nb,style='C'), method="eigen",
                           quiet=FALSE); summary(spat.err.immig)
## 
## Spatial autoregressive error model
## 
## Jacobian calculated using neighbourhood matrix eigenvalues
## Computing eigenvalues ...
## 
## lambda: -0.6478654  function: -51.85541  Jacobian: -2.10781  SSE: 22.33468 
## lambda: -0.08548168  function: -47.84376  Jacobian: -0.0385715  SSE: 20.59826 
## lambda: 0.2620906  function: -47.83341  Jacobian: -0.4034806  SSE: 20.2787 
## lambda: 0.09021095  function: -47.58365  Jacobian: -0.04491613  SSE: 20.37085 
## lambda: 0.09006783  function: -47.58366  Jacobian: -0.04477178  SSE: 20.37097 
## lambda: 0.09335054  function: -47.58356  Jacobian: -0.04814291  SSE: 20.36803 
## lambda: 0.1578035  function: -47.61883  Jacobian: -0.1404927  SSE: 20.31965 
## lambda: 0.0934995  function: -47.58356  Jacobian: -0.04829886  SSE: 20.3679 
## lambda: 0.09349818  function: -47.58356  Jacobian: -0.04829749  SSE: 20.3679 
## lambda: 0.09349813  function: -47.58356  Jacobian: -0.04829743  SSE: 20.3679 
## lambda: 0.09349816  function: -47.58356  Jacobian: -0.04829747  SSE: 20.3679 
## lambda: 0.09349868  function: -47.58356  Jacobian: -0.04829801  SSE: 20.3679 
## lambda: 0.09349837  function: -47.58356  Jacobian: -0.04829769  SSE: 20.3679 
## lambda: 0.09349826  function: -47.58356  Jacobian: -0.04829756  SSE: 20.3679 
## lambda: 0.09349821  function: -47.58356  Jacobian: -0.04829752  SSE: 20.3679 
## lambda: 0.09349819  function: -47.58356  Jacobian: -0.0482975  SSE: 20.3679 
## lambda: 0.09349817  function: -47.58356  Jacobian: -0.04829748  SSE: 20.3679 
## lambda: 0.09349818  function: -47.58356  Jacobian: -0.04829749  SSE: 20.3679
## 
## Call:errorsarlm(formula = immig0511 ~ pubIdeolCCES, data = imm.data, 
##     listw = nb2listw(state.nb, style = "C"), method = "eigen", 
##     quiet = FALSE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.33868 -0.51886  0.02905  0.47115  1.37927 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value        Pr(>|z|)
## (Intercept)  0.922768   0.170728  5.4049 0.0000000648425
## pubIdeolCCES 0.065708   0.010713  6.1335 0.0000000008597
## 
## Lambda: 0.093498, LR test value: 0.14407, p-value: 0.70427
## Asymptotic standard error: 0.2068
##     z-value: 0.45213, p-value: 0.65118
## Wald statistic: 0.20442, p-value: 0.65118
## 
## Log likelihood: -47.58356 for error model
## ML residual variance (sigma squared): 0.42433, (sigma: 0.65141)
## Number of observations: 48 
## Number of parameters estimated: 4 
## AIC: 103.17, (AIC for lm: 101.31)
#Fit spatial lag model
spat.lag.immig<-lagsarlm(immig0511~pubIdeolCCES, data=imm.data,
                         nb2listw(state.nb,style='C'), method="eigen",
                         quiet=FALSE); summary(spat.lag.immig)
## 
## Spatial lag model
## Jacobian calculated using neighbourhood matrix eigenvalues
## Computing eigenvalues ...
## 
## rho:  -0.6478654     function value:  -57.83021 
## rho:  -0.08548168    function value:  -48.34439 
## rho:  0.2620906  function value:  -47.17023 
## rho:  0.2022487  function value:  -47.04967 
## rho:  0.1784245  function value:  -47.04129 
## rho:  0.1814846  function value:  -47.04113 
## rho:  0.1813312  function value:  -47.04112 
## rho:  0.1813374  function value:  -47.04112 
## rho:  0.1813374  function value:  -47.04112 
## rho:  0.1813374  function value:  -47.04112 
## rho:  0.1813374  function value:  -47.04112 
## rho:  0.1813374  function value:  -47.04112
## 
## Call:lagsarlm(formula = immig0511 ~ pubIdeolCCES, data = imm.data, 
##     listw = nb2listw(state.nb, style = "C"), method = "eigen", 
##     quiet = FALSE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.34985 -0.41520  0.08061  0.49837  1.40920 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value      Pr(>|z|)
## (Intercept)  0.862126   0.171690  5.0214 0.00000051292
## pubIdeolCCES 0.061133   0.011192  5.4621 0.00000004707
## 
## Rho: 0.18134, LR test value: 1.2289, p-value: 0.26761
## Asymptotic standard error: 0.15927
##     z-value: 1.1386, p-value: 0.25489
## Wald statistic: 1.2963, p-value: 0.25489
## 
## Log likelihood: -47.04112 for lag model
## ML residual variance (sigma squared): 0.41246, (sigma: 0.64223)
## Number of observations: 48 
## Number of parameters estimated: 4 
## AIC: 102.08, (AIC for lm: 101.31)
## LM test for residual autocorrelation
## test value: 0.54053, p-value: 0.46221

0.0.4 Point-Referenced Data Models

###Replication of Scallop example in R from chapter 2 of Banerjee, Carlin, & Gelfand 2004###

rm(list=ls())

#input data
myscallops <- read.table("myscallops.txt", header=T)

#load libraries
library(akima) #for drawing image/surface plots
library(geoR) #for variogram analyses
## --------------------------------------------------------------
##  Analysis of Geostatistical Data
##  For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
##  geoR version 1.7-5.2.1 (built on 2016-05-02) is now loaded
## --------------------------------------------------------------
library(scatterplot3d) #for drop line plot

#create an image plot with contours
int.scp <- akima::interp(x=myscallops$long, 
                         y=myscallops$lat, 
                         z=myscallops$lgcatch)
image(int.scp, xlim=range(myscallops$long), 
      ylim=range(myscallops$lat))
contour(int.scp, add=T)

#create a surface plot
persp(int.scp, xlim=range(myscallops$long), 
      ylim=range(myscallops$lat))

#drop line plot
scatterplot3d(x=myscallops$long, y=myscallops$lat, z=myscallops$lgcatch,
              pch=20, highlight.3d=TRUE, type='h')

#define a "geodata" object
obj <- cbind(myscallops$long, myscallops$lat, myscallops$lgcatch)
scallops.geo <- geoR::as.geodata(obj, coords.col=1:2, data.col=3)

#create an empirical semivariogram
scallops.var <- geoR::variog(scallops.geo, estimator.type="classical")
## variog: computing omnidirectional variogram
scallops.var
## $u
##  [1] 0.1115405 0.3346215 0.5577025 0.7807835 1.0038645 1.2269455 1.4500265
##  [8] 1.6731075 1.8961885 2.1192695 2.3423505 2.5654315 2.7885125
## 
## $v
##  [1] 3.048744 4.361375 4.520836 4.924096 5.697811 5.175957 4.888780
##  [8] 4.345750 4.455269 4.224356 4.460676 5.338326 7.391131
## 
## $n
##  [1]  669 1538 1566 1534 1410 1267 1011  724  465  329  218  109   37
## 
## $sd
##  [1] 5.736655 6.292244 5.600070 6.288254 7.001918 6.329819 6.570810
##  [8] 5.975687 5.815976 4.899078 5.451433 5.928083 7.090025
## 
## $bins.lim
##  [1] 0.000000000001 0.223081002026 0.446162004053 0.669243006079
##  [5] 0.892324008105 1.115405010132 1.338486012158 1.561567014185
##  [9] 1.784648016211 2.007729018237 2.230810020264 2.453891022290
## [13] 2.676972024316 2.900053026343
## 
## $ind.bin
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 
## $var.mark
## [1] 4.720595
## 
## $beta.ols
## [1] 3.482623
## 
## $output.type
## [1] "bin"
## 
## $max.dist
## [1] 2.900053
## 
## $estimator.type
## [1] "classical"
## 
## $n.data
## [1] 148
## 
## $lambda
## [1] 1
## 
## $trend
## [1] "cte"
## 
## $pairs.min
## [1] 2
## 
## $nugget.tolerance
## [1] 1e-12
## 
## $direction
## [1] "omnidirectional"
## 
## $tolerance
## [1] "none"
## 
## $uvec
##  [1] 0.1115405 0.3346215 0.5577025 0.7807835 1.0038645 1.2269455 1.4500265
##  [8] 1.6731075 1.8961885 2.1192695 2.3423505 2.5654315 2.7885125
## 
## $call
## geoR::variog(geodata = scallops.geo, estimator.type = "classical")
## 
## attr(,"class")
## [1] "variogram"
#create a robust empirical semivariogram
scallops.var.robust <- geoR::variog(scallops.geo, estimator.type="modulus")
## variog: computing omnidirectional variogram
#plot the two forms of the semivariogram
par(mfrow=c(2,1)) #align separate plots on the same space, two rows by one column
plot(scallops.var)
plot(scallops.var.robust)

#Fit a parametric variogram model (assume an exponential model) using weighted least squares
scallops.var.fit <- variofit(scallops.var.robust, 
                             ini.cov.pars = c(6.0,1.0), 
                             cov.model="exponential", 
                             fix.nugget=FALSE, nugget=1.0)
## variofit: covariance model used is exponential 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim
scallops.var.fit
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
##   tausq sigmasq     phi 
##  0.0000  5.1110  0.2065 
## Practical Range with cor=0.05 for asymptotic range: 0.6184923
## 
## variofit: minimised weighted sum of squares = 4017.547
#Fit the exponential model using MLE
scallops.lik.fit <- likfit(scallops.geo, ini.cov.pars=c(6.0,1.0), 
                           cov.model = "exponential", trend = "cte",
                           fix.nugget = FALSE, nugget = 1.0, 
                           nospatial = TRUE, lik.method = "ML") 
## kappa not used for the exponential correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.
#note: "lik.method" is a correction of a typo from the text
scallops.lik.fit
## likfit: estimated model parameters:
##     beta    tausq  sigmasq      phi 
## "2.3748" "0.0947" "5.7675" "0.2338" 
## Practical Range with cor=0.05 for asymptotic range: 0.700551
## 
## likfit: maximised log-likelihood = -285.9
#CLASSICAL KRIGING
#First we must create locations to predict#
#Source: http://sudipto.bol.ucla.edu/Software/Book.R
u.loci <- rbind(c(-73.0,39.5),c(-72.5,40.0),c(-72,39.0))

#visualize where we're predicting
plot(y=myscallops$lat, x=myscallops$long, pch='.')
points(u.loci, pch='+', col='red')

scallops.classical <- geoR::krige.conv(geodata = scallops.geo, 
                                       locations = u.loci, borders = NULL,
                                       krige = krige.control(cov.pars=c(5.7675,0.2338))) 
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
#cov.pars=c(sigmasq,phi)
scallops.classical$predict
## [1] 5.675546 4.972033 2.067120
#Bayesian Kriging Model--Estimation and Prediction
scallops.bayes2 <- geoR::krige.bayes(scallops.geo,
    locations = u.loci, #defaults to locations="no" 
    borders = NULL,
    model = model.control(trend.d = "cte",
        cov.model = "exponential"),
    prior = prior.control(beta.prior = "flat",
        sigmasq.prior = "reciprocal",
        tausq.rel.prior = "uniform",
        tausq.rel.discrete = seq(from=0.0, to=1.0, by=0.01)))
## Warning in geoR::krige.bayes(scallops.geo, locations = u.loci, borders =
## NULL, : .Random.seed not initialised. Creating it with runif(1)
## krige.bayes: model with constant mean
## krige.bayes: computing the discrete posterior of phi/tausq.rel
## krige.bayes: argument `phi.discrete` not provided, using default values
## krige.bayes: computing the posterior probabilities.
##              Number of parameter sets:  5050 
## 1, 101, 201, 301, 401, 501, 601, 701, 801, 901, 1001, 1101, 1201, 1301, 1401, 1501, 1601, 1701, 1801, 1901, 2001, 2101, 2201, 2301, 2401, 2501, 2601, 2701, 2801, 2901, 3001, 3101, 3201, 3301, 3401, 3501, 3601, 3701, 3801, 3901, 4001, 4101, 4201, 4301, 4401, 4501, 4601, 4701, 4801, 4901, 5001, 
## 
## krige.bayes: sampling from posterior distribution
## krige.bayes: sample from the (joint) posterior of phi and tausq.rel
##                 [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
## phi        0.2320042 0.3480064 0.4640085 0.5800106 0.6960127 0.8120148
## tausq.rel  0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## frequency 18.0000000 7.0000000 5.0000000 8.0000000 9.0000000 4.0000000
##               [,7]     [,8]     [,9]    [,10]    [,11]   [,12]    [,13]
## phi       0.928017 1.044019 1.160021 1.276023 1.392025 1.62403 1.740032
## tausq.rel 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000
## frequency 4.000000 3.000000 3.000000 4.000000 4.000000 3.00000 2.000000
##              [,14]    [,15]    [,16]   [,17]    [,18]    [,19]    [,20]
## phi       1.856034 1.972036 2.088038 2.20404 2.320042 2.436045 2.552047
## tausq.rel 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000
## frequency 2.000000 4.000000 2.000000 3.00000 5.000000 1.000000 4.000000
##              [,21]    [,22]    [,23]    [,24]    [,25]    [,26]    [,27]
## phi       2.668049 2.784051 2.900053 3.016055 3.132057 3.364062 3.480064
## tausq.rel 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
## frequency 4.000000 2.000000 1.000000 3.000000 1.000000 1.000000 1.000000
##              [,28]    [,29]   [,30]    [,31]    [,32]    [,33]    [,34]
## phi       3.596066 3.712068 3.82807 3.944072 4.060074 4.176076 4.292078
## tausq.rel 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000 0.000000
## frequency 2.000000 2.000000 4.00000 1.000000 5.000000 1.000000 2.000000
##              [,35]    [,36]    [,37]    [,38]    [,39]    [,40]    [,41]
## phi       4.408081 4.524083 4.640085 4.872089 4.988091 5.104093 5.220095
## tausq.rel 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
## frequency 1.000000 2.000000 3.000000 2.000000 2.000000 2.000000 1.000000
##              [,42]  [,43]    [,44]    [,45]    [,46]      [,47]      [,48]
## phi       5.336098 5.4521 5.568102 5.684104 5.800106  0.2320042  0.3480064
## tausq.rel 0.000000 0.0000 0.000000 0.000000 0.000000  0.0100000  0.0100000
## frequency 2.000000 1.0000 1.000000 4.000000 1.000000 12.0000000 19.0000000
##               [,49]     [,50]     [,51]     [,52]    [,53]    [,54]
## phi       0.4640085 0.5800106 0.6960127 0.8120148 0.928017 1.044019
## tausq.rel 0.0100000 0.0100000 0.0100000 0.0100000 0.010000 0.010000
## frequency 5.0000000 9.0000000 8.0000000 3.0000000 6.000000 7.000000
##              [,55]    [,56]    [,57]    [,58]   [,59]    [,60]    [,61]
## phi       1.160021 1.276023 1.392025 1.508028 1.62403 1.740032 1.856034
## tausq.rel 0.010000 0.010000 0.010000 0.010000 0.01000 0.010000 0.010000
## frequency 2.000000 2.000000 2.000000 2.000000 7.00000 5.000000 7.000000
##              [,62]    [,63]   [,64]    [,65]    [,66]    [,67]    [,68]
## phi       1.972036 2.088038 2.20404 2.320042 2.436045 2.552047 2.668049
## tausq.rel 0.010000 0.010000 0.01000 0.010000 0.010000 0.010000 0.010000
## frequency 5.000000 6.000000 6.00000 4.000000 3.000000 4.000000 1.000000
##              [,69]    [,70]    [,71]    [,72]    [,73]    [,74]    [,75]
## phi       2.784051 2.900053 3.016055 3.132057 3.248059 3.364062 3.480064
## tausq.rel 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000
## frequency 4.000000 6.000000 1.000000 1.000000 3.000000 2.000000 5.000000
##              [,76]    [,77]   [,78]    [,79]    [,80]    [,81]    [,82]
## phi       3.596066 3.712068 3.82807 4.060074 4.176076 4.292078 4.408081
## tausq.rel 0.010000 0.010000 0.01000 0.010000 0.010000 0.010000 0.010000
## frequency 4.000000 2.000000 2.00000 2.000000 4.000000 2.000000 5.000000
##              [,83]    [,84]    [,85]    [,86]    [,87]    [,88]  [,89]
## phi       4.524083 4.640085 4.872089 4.988091 5.220095 5.336098 5.4521
## tausq.rel 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.0100
## frequency 2.000000 3.000000 1.000000 1.000000 1.000000 5.000000 1.0000
##              [,90]    [,91]    [,92]      [,93]      [,94]      [,95]
## phi       5.568102 5.684104 5.800106  0.2320042  0.3480064  0.4640085
## tausq.rel 0.010000 0.010000 0.010000  0.0200000  0.0200000  0.0200000
## frequency 3.000000 1.000000 1.000000 16.0000000 11.0000000 16.0000000
##                [,96]      [,97]     [,98]     [,99]   [,100]   [,101]
## phi        0.5800106  0.6960127 0.8120148  0.928017 1.044019 1.160021
## tausq.rel  0.0200000  0.0200000 0.0200000  0.020000 0.020000 0.020000
## frequency 10.0000000 11.0000000 6.0000000 12.000000 3.000000 3.000000
##             [,102]   [,103]   [,104]  [,105]   [,106]   [,107]   [,108]
## phi       1.276023 1.392025 1.508028 1.62403 1.740032 1.856034 1.972036
## tausq.rel 0.020000 0.020000 0.020000 0.02000 0.020000 0.020000 0.020000
## frequency 6.000000 5.000000 2.000000 2.00000 8.000000 5.000000 5.000000
##             [,109]  [,110]   [,111]   [,112]   [,113]   [,114]   [,115]
## phi       2.088038 2.20404 2.320042 2.436045 2.552047 2.668049 2.900053
## tausq.rel 0.020000 0.02000 0.020000 0.020000 0.020000 0.020000 0.020000
## frequency 2.000000 3.00000 2.000000 5.000000 2.000000 2.000000 2.000000
##             [,116]   [,117]   [,118]   [,119]    [,120]     [,121]
## phi       3.016055 3.132057 3.596066 4.640085 0.1160021  0.2320042
## tausq.rel 0.020000 0.020000 0.020000 0.020000 0.0300000  0.0300000
## frequency 2.000000 1.000000 1.000000 1.000000 2.0000000 18.0000000
##               [,122]     [,123]    [,124]    [,125]    [,126]   [,127]
## phi        0.3480064  0.4640085 0.5800106 0.6960127 0.8120148 0.928017
## tausq.rel  0.0300000  0.0300000 0.0300000 0.0300000 0.0300000 0.030000
## frequency 17.0000000 14.0000000 5.0000000 7.0000000 4.0000000 4.000000
##             [,128]   [,129]   [,130]   [,131]   [,132]  [,133]   [,134]
## phi       1.044019 1.160021 1.276023 1.392025 1.508028 1.62403 1.740032
## tausq.rel 0.030000 0.030000 0.030000 0.030000 0.030000 0.03000 0.030000
## frequency 3.000000 6.000000 2.000000 4.000000 1.000000 4.00000 3.000000
##             [,135]   [,136]  [,137]   [,138]   [,139]   [,140]   [,141]
## phi       1.856034 1.972036 2.20404 2.552047 2.784051 3.364062 4.872089
## tausq.rel 0.030000 0.030000 0.03000 0.030000 0.030000 0.030000 0.030000
## frequency 1.000000 1.000000 1.00000 1.000000 2.000000 1.000000 1.000000
##               [,142]     [,143]     [,144]    [,145]    [,146]    [,147]
## phi        0.2320042  0.3480064  0.4640085 0.5800106 0.6960127 0.8120148
## tausq.rel  0.0400000  0.0400000  0.0400000 0.0400000 0.0400000 0.0400000
## frequency 13.0000000 11.0000000 12.0000000 3.0000000 5.0000000 7.0000000
##             [,148]   [,149]   [,150]   [,151]   [,152]  [,153]   [,154]
## phi       0.928017 1.044019 1.160021 1.276023 1.392025 1.62403 1.972036
## tausq.rel 0.040000 0.040000 0.040000 0.040000 0.040000 0.04000 0.040000
## frequency 6.000000 4.000000 1.000000 4.000000 1.000000 2.00000 1.000000
##             [,155]   [,156]  [,157]    [,158]     [,159]    [,160]
## phi       2.436045 2.784051 3.82807 0.2320042  0.3480064 0.4640085
## tausq.rel 0.040000 0.040000 0.04000 0.0500000  0.0500000 0.0500000
## frequency 1.000000 1.000000 1.00000 9.0000000 14.0000000 7.0000000
##              [,161]    [,162]    [,163]   [,164]   [,165]   [,166]
## phi       0.5800106 0.6960127 0.8120148 0.928017 1.044019 1.160021
## tausq.rel 0.0500000 0.0500000 0.0500000 0.050000 0.050000 0.050000
## frequency 4.0000000 9.0000000 3.0000000 3.000000 3.000000 4.000000
##             [,167]   [,168]  [,169]   [,170]     [,171]    [,172]
## phi       1.276023 1.392025 1.62403 1.740032  0.2320042 0.3480064
## tausq.rel 0.050000 0.050000 0.05000 0.050000  0.0600000 0.0600000
## frequency 1.000000 2.000000 1.00000 1.000000 14.0000000 9.0000000
##               [,173]    [,174]    [,175]   [,176]   [,177]   [,178]
## phi        0.4640085 0.5800106 0.8120148 0.928017 1.044019 1.276023
## tausq.rel  0.0600000 0.0600000 0.0600000 0.060000 0.060000 0.060000
## frequency 10.0000000 8.0000000 5.0000000 3.000000 4.000000 2.000000
##              [,179]    [,180]     [,181]    [,182]    [,183]    [,184]
## phi       0.1160021 0.2320042  0.3480064 0.4640085 0.5800106 0.6960127
## tausq.rel 0.0700000 0.0700000  0.0700000 0.0700000 0.0700000 0.0700000
## frequency 1.0000000 5.0000000 18.0000000 3.0000000 5.0000000 3.0000000
##             [,185]    [,186]    [,187]    [,188]    [,189]    [,190]
## phi       1.972036 0.2320042 0.3480064 0.4640085 0.5800106 0.6960127
## tausq.rel 0.070000 0.0800000 0.0800000 0.0800000 0.0800000 0.0800000
## frequency 1.000000 7.0000000 9.0000000 7.0000000 4.0000000 2.0000000
##              [,191]   [,192]   [,193]   [,194]   [,195]    [,196]
## phi       0.8120148 0.928017 1.044019 1.276023 1.740032 0.2320042
## tausq.rel 0.0800000 0.080000 0.080000 0.080000 0.080000 0.0900000
## frequency 1.0000000 3.000000 1.000000 1.000000 1.000000 9.0000000
##              [,197]    [,198]    [,199]    [,200]    [,201]   [,202]
## phi       0.3480064 0.4640085 0.5800106 0.6960127 0.8120148 1.276023
## tausq.rel 0.0900000 0.0900000 0.0900000 0.0900000 0.0900000 0.090000
## frequency 8.0000000 7.0000000 6.0000000 2.0000000 2.0000000 1.000000
##              [,203]    [,204]    [,205]    [,206]    [,207]   [,208]
## phi       0.2320042 0.3480064 0.4640085 0.5800106 0.8120148 0.928017
## tausq.rel 0.1000000 0.1000000 0.1000000 0.1000000 0.1000000 0.100000
## frequency 2.0000000 6.0000000 6.0000000 2.0000000 1.0000000 1.000000
##              [,209]    [,210]    [,211]    [,212]    [,213]    [,214]
## phi       0.1160021 0.2320042 0.3480064 0.4640085 0.5800106 0.8120148
## tausq.rel 0.1100000 0.1100000 0.1100000 0.1100000 0.1100000 0.1100000
## frequency 1.0000000 3.0000000 4.0000000 2.0000000 3.0000000 1.0000000
##              [,215]    [,216]    [,217]    [,218]    [,219]    [,220]
## phi       0.2320042 0.3480064 0.4640085 0.5800106 0.6960127 0.8120148
## tausq.rel 0.1200000 0.1200000 0.1200000 0.1200000 0.1200000 0.1200000
## frequency 7.0000000 5.0000000 1.0000000 4.0000000 1.0000000 1.0000000
##              [,221]    [,222]    [,223]    [,224]    [,225]    [,226]
## phi       0.2320042 0.3480064 0.4640085 0.5800106 0.2320042 0.3480064
## tausq.rel 0.1300000 0.1300000 0.1300000 0.1300000 0.1400000 0.1400000
## frequency 2.0000000 6.0000000 3.0000000 2.0000000 3.0000000 4.0000000
##              [,227]    [,228]    [,229]    [,230]    [,231]    [,232]
## phi       0.4640085 0.5800106 0.8120148 0.2320042 0.3480064 0.4640085
## tausq.rel 0.1400000 0.1400000 0.1400000 0.1500000 0.1500000 0.1500000
## frequency 1.0000000 1.0000000 1.0000000 2.0000000 4.0000000 1.0000000
##              [,233]    [,234]    [,235]    [,236]    [,237]    [,238]
## phi       0.2320042 0.3480064 0.6960127 0.2320042 0.3480064 0.4640085
## tausq.rel 0.1600000 0.1600000 0.1600000 0.1700000 0.1700000 0.1700000
## frequency 2.0000000 4.0000000 1.0000000 2.0000000 1.0000000 1.0000000
##              [,239]    [,240]    [,241]    [,242]    [,243]    [,244]
## phi       0.2320042 0.4640085 0.2320042 0.3480064 0.2320042 0.2320042
## tausq.rel 0.1800000 0.1800000 0.1900000 0.1900000 0.2000000 0.2100000
## frequency 6.0000000 4.0000000 2.0000000 4.0000000 1.0000000 1.0000000
##              [,245]    [,246]    [,247]    [,248]    [,249]    [,250]
## phi       0.5800106 0.2320042 0.2320042 0.4640085 0.3480064 0.4640085
## tausq.rel 0.2200000 0.2300000 0.2500000 0.2600000 0.2700000 0.2700000
## frequency 1.0000000 1.0000000 1.0000000 1.0000000 3.0000000 1.0000000
##              [,251]    [,252]    [,253]    [,254]    [,255]    [,256]
## phi       0.5800106 0.4640085 0.2320042 0.4640085 0.2320042 0.2320042
## tausq.rel 0.2700000 0.2900000 0.3100000 0.3600000 0.4100000 0.5300000
## frequency 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## 
## krige.bayes: starting prediction at the provided locations
## krige.bayes: phi/tausq.rel samples for the predictive are same as for the posterior 
## krige.bayes: computing moments of the predictive distribution
## krige.bayes: sampling from the predictive
##              Number of parameter sets:  256 
## 1, 11, 21, 31, 41, 51, 61, 71, 81, 91, 101, 111, 121, 131, 141, 151, 161, 171, 181, 191, 201, 211, 221, 231, 241, 251, 
## krige.bayes: preparing summaries of the predictive distribution
# Projected number of ads at three sites from Bayesian model
#Note: sites are set by: u.loci <- rbind(c(-73.0,39.5),c(-72.5,40.0),c(-72,39.0))

scallops.bayes2$predictive$mean
## [1]  5.6620856  5.0163805 -0.3789692
#Form the quantiles to interpret the model results
out <- scallops.bayes2$posterior
out <- out$sample
beta.qnt <- quantile(out$beta, c(0.50,0.025,0.975))
phi.qnt <-quantile(out$phi, c(0.50,0.025,0.975))
sigmasq.qnt <- quantile(out$sigmasq, c(0.50,0.025,0.975))
tausq.rel.qnt <- quantile(out$tausq.rel, c(0.50, 0.025, 0.975))
beta.qnt
##       50%      2.5%     97.5% 
##  1.821189 -6.790915  7.366864
phi.qnt
##       50%      2.5%     97.5% 
## 0.5800106 0.2320042 4.9880912
sigmasq.qnt
##       50%      2.5%     97.5% 
## 10.994736  4.226533 94.983472
tausq.rel.qnt
##   50%  2.5% 97.5% 
##  0.03  0.00  0.18
######################################################################

#R code for mapping locations of the 1974 pilot data
#clean up
rm(list=ls())

#load libraries
library(spatstat)
## Loading required package: spatstat.data
## Loading required package: nlme
## Loading required package: rpart
## 
## spatstat 1.57-1       (nickname: 'Cartoon Physics') 
## For an introduction to spatstat, type 'beginner'
library(maptools)
library(rgdal)
library(maps)
library(foreign)
library(mgcv)
## This is mgcv 1.8-26. For overview type 'help("mgcv-package")'.
library(rgeos)
## rgeos version: 0.4-2, (SVN revision 581)
##  GEOS runtime version: 3.6.1-CAPI-1.10.1 
##  Linking to sp version: 1.3-1 
##  Polygon checking: TRUE
library(spdep)
library(ggmap)
## Loading required package: ggplot2
checkCRSArgs <- function(uprojargs) {
  .Call("checkCRSArgs", uprojargs, PACKAGE="rgdal")
  }

#options
options(scipen=8)

#load the 1974 pilot data
pilot<-read.csv("pilot74.csv")

#visualize the locations
us.states<-map("state",fill=FALSE, plot=TRUE, resolution=0)

IDs <- us.states$names
us.states.1<-map(database="state",region=IDs, fill=TRUE, 
                 plot=TRUE, interior=TRUE, exact=TRUE, resolution=0)
a.1<-maptools::map2SpatialPolygons(us.states.1, IDs=IDs, 
                                   proj4string=CRS("+proj=longlat"))
us.map.1<-spTransform(a.1,CRS("+init=esri:102003"))
###USEFUL WEBSITES###
#  http://spatialreference.org/ref/esri/102003/
#  https://www.nhgis.org/

plot(us.map.1,xlim=us.map.1@bbox[1,],ylim=us.map.1@bbox[2,])
axis(1);axis(2);box()
points(y=pilot$northings,x=pilot$eastings,col='red',pch='.')

0.0.5 Basics of Areal Data Models

#Estimating CAR models with world democracy data
#clean up
rm(list=ls())

#preferences
options(scipen=8)

#load required libraries
library(spdep) 
library(maps)
library(foreign)
library(rgdal)

#Load Data
# sldv is the data.frame
# mdd2 is the minimum distance data.frame
source("https://spia.uga.edu/faculty_pages/monogan/teaching/spatial/chapter1data.R")
#fix(sldv)
#names(sldv)
#fix(mdd2)
#names(mdd2)

# Create a neighbor minimum distance matrix
# from the dyadic minimum distance data
# (Note that the cshapes library provides a 
# convienent way to create a matrix directly)
mddmat <- matrix(9999,ncol=dim(sldv)[1],nrow=dim(sldv)[1]) 
#We choose 9999km as the default because it's a big number that indicates "not a neighbor."
dimnames(mddmat) <- list(c(sldv$tla),c(sldv$tla))
for(i in 1:dim(mdd2)[1]){
    mddmat[mdd2$ida[i],mdd2$idb[i]] <-  mdd2$mindist[i] 
    #super cool stuff: call text entries to set matrix cells
}
mddmat[1:12,1:12]
##      AFG  ALB  ALG  ANG  ARG  ARM  AUL  AUS  AZE  BAH  BEL  BEN
## AFG 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999
## ALB 9999 9999 9999 9999 9999 9999 9999  554 9999 9999 9999 9999
## ALG 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999  763
## ANG 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999
## ARG 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999
## ARM 9999 9999 9999 9999 9999 9999 9999 9999    0 9999 9999 9999
## AUL 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999
## AUS 9999  554 9999 9999 9999 9999 9999 9999 9999 9999  369 9999
## AZE 9999 9999 9999 9999 9999    0 9999 9999 9999 9999 9999 9999
## BAH 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999
## BEL 9999 9999 9999 9999 9999 9999 9999  369 9999 9999 9999 9999
## BEN 9999 9999  763 9999 9999 9999 9999 9999 9999 9999 9999 9999
# Create a binary matrix with a 200 km threshold
m200mat <- mddmat
m200mat[m200mat<=200] <- 1
m200mat[m200mat>200] <- 0
m200mat[1:12,1:12]
##     AFG ALB ALG ANG ARG ARM AUL AUS AZE BAH BEL BEN
## AFG   0   0   0   0   0   0   0   0   0   0   0   0
## ALB   0   0   0   0   0   0   0   0   0   0   0   0
## ALG   0   0   0   0   0   0   0   0   0   0   0   0
## ANG   0   0   0   0   0   0   0   0   0   0   0   0
## ARG   0   0   0   0   0   0   0   0   0   0   0   0
## ARM   0   0   0   0   0   0   0   0   1   0   0   0
## AUL   0   0   0   0   0   0   0   0   0   0   0   0
## AUS   0   0   0   0   0   0   0   0   0   0   0   0
## AZE   0   0   0   0   0   1   0   0   0   0   0   0
## BAH   0   0   0   0   0   0   0   0   0   0   0   0
## BEL   0   0   0   0   0   0   0   0   0   0   0   0
## BEN   0   0   0   0   0   0   0   0   0   0   0   0
# Create a row standardized weights matrix
wmat <- m200mat/apply(m200mat,1,sum)
wmat[1:12,1:12]
##     AFG ALB ALG ANG ARG       ARM AUL AUS AZE BAH BEL BEN
## AFG   0   0   0   0   0 0.0000000   0   0 0.0   0   0   0
## ALB   0   0   0   0   0 0.0000000   0   0 0.0   0   0   0
## ALG   0   0   0   0   0 0.0000000   0   0 0.0   0   0   0
## ANG   0   0   0   0   0 0.0000000   0   0 0.0   0   0   0
## ARG   0   0   0   0   0 0.0000000   0   0 0.0   0   0   0
## ARM   0   0   0   0   0 0.0000000   0   0 0.2   0   0   0
## AUL   0   0   0   0   0 0.0000000   0   0 0.0   0   0   0
## AUS   0   0   0   0   0 0.0000000   0   0 0.0   0   0   0
## AZE   0   0   0   0   0 0.1428571   0   0 0.0   0   0   0
## BAH   0   0   0   0   0 0.0000000   0   0 0.0   0   0   0
## BEL   0   0   0   0   0 0.0000000   0   0 0.0   0   0   0
## BEN   0   0   0   0   0 0.0000000   0   0 0.0   0   0   0
# calculate the spatial lag of democracy
democracy.spatial.lag <- as.vector(wmat%*%sldv$democracy)
democracy.spatial.lag[1:12]
##  [1]  -4.5714286   7.0000000  -0.4285714   2.2000000   8.6000000
##  [6]   3.0000000   8.5000000   8.5454545   1.7142857 -10.0000000
## [11]   9.8000000   2.4000000
# Create a weights list object (a list object is required for spdep)
cmat.lw <- mat2listw(m200mat, row.names=row.names(m200mat))

# Extract a neighbor list
cmat.nb <- cmat.lw$neighbours

#Fit Conditional Autoregressive Error Model
sldv.car<-spautolm(democracy ~ log(gdp.2002/population), data=sldv, listw=nb2listw(cmat.nb), family='CAR'); summary(sldv.car)
## Warning in spautolm(democracy ~ log(gdp.2002/population), data = sldv,
## listw = nb2listw(cmat.nb), : Non-symmetric spatial weights in CAR model
## 
## Call: 
## spautolm(formula = democracy ~ log(gdp.2002/population), data = sldv, 
##     listw = nb2listw(cmat.nb), family = "CAR")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -16.84735  -3.37506   0.44561   2.94816  12.81756 
## 
## Coefficients: 
##                          Estimate Std. Error z value Pr(>|z|)
## (Intercept)              -5.17337    3.32570 -1.5556 0.119810
## log(gdp.2002/population)  1.19528    0.38461  3.1078 0.001885
## 
## Lambda: 0.91466 LR test value: 50.195 p-value: 0.0000000000013922 
## Numerical Hessian standard error of lambda: 0.061333 
## 
## Log likelihood: -488.5194 
## ML residual variance (sigma squared): 24.496, (sigma: 4.9493)
## Number of observations: 158 
## Number of parameters estimated: 4 
## AIC: 985.04
#Fit spatial error model from Ward & Gleditsch
sldv.err<-errorsarlm(democracy ~ log(gdp.2002/population), data=sldv, listw=nb2listw(cmat.nb), method='eigen'); summary(sldv.err)
## 
## Call:
## errorsarlm(formula = democracy ~ log(gdp.2002/population), data = sldv, 
##     listw = nb2listw(cmat.nb), method = "eigen")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.9490  -4.2318   1.4473   3.4322  11.5587 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                          Estimate Std. Error z value  Pr(>|z|)
## (Intercept)              -7.48651    3.07263 -2.4365 0.0148296
## log(gdp.2002/population)  1.38711    0.37911  3.6589 0.0002533
## 
## Lambda: 0.58188, LR test value: 44.18, p-value: 0.000000000029953
## Asymptotic standard error: 0.076529
##     z-value: 7.6034, p-value: 2.8866e-14
## Wald statistic: 57.811, p-value: 2.8866e-14
## 
## Log likelihood: -491.5267 for error model
## ML residual variance (sigma squared): 27.14, (sigma: 5.2096)
## Number of observations: 158 
## Number of parameters estimated: 4 
## AIC: 991.05, (AIC for lm: 1033.2)
#####################################################################
#Areal models and diagnostic tests with immigrant policy data

#Load Data
imm.data.0<-read.dta("stateImmig0511.dta")
#subset
imm.data<-subset(imm.data.0,State!="ALASKA" & State!="HAWAII")

#call data from "maps", create neighbor matrix
us.states<-map("state", fill=TRUE, plot=TRUE, resolution=0)

IDs <- sapply(strsplit(us.states$names, ":"), function(x) x[1])
us.map.0<-map2SpatialPolygons(us.states, IDs=IDs)

#eliminate Washington, DC
us.map<-us.map.0[-8]

#create neighbor list from state map
state.nb<-poly2nb(us.map)

#Basic linear model
mod.immigrant<-lm(immig0511~pubIdeolCCES+sqrt(squireProfess)+repUnif+
                    demUnif+pcgsp1000+termLimits+changeForeign, data=imm.data); summary(mod.immigrant)
## 
## Call:
## lm(formula = immig0511 ~ pubIdeolCCES + sqrt(squireProfess) + 
##     repUnif + demUnif + pcgsp1000 + termLimits + changeForeign, 
##     data = imm.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.3148 -0.3415  0.1175  0.4027  0.7816 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)   
## (Intercept)         -1.217430   0.725108  -1.679  0.10096   
## pubIdeolCCES         0.037442   0.012468   3.003  0.00459 **
## sqrt(squireProfess)  1.643635   0.766027   2.146  0.03802 * 
## repUnif              0.038726   0.037425   1.035  0.30700   
## demUnif              0.021630   0.039944   0.542  0.59116   
## pcgsp1000            0.032330   0.011810   2.737  0.00919 **
## termLimits          -0.244108   0.189977  -1.285  0.20621   
## changeForeign       -0.004057   0.001470  -2.760  0.00868 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5746 on 40 degrees of freedom
## Multiple R-squared:  0.6572, Adjusted R-squared:  0.5972 
## F-statistic: 10.95 on 7 and 40 DF,  p-value: 0.0000001241
#tests on this linear model
tests<-lm.LMtests(mod.immigrant, 
                  listw=nb2listw(state.nb, style='W'),test="all");tests
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = immig0511 ~ pubIdeolCCES + sqrt(squireProfess)
## + repUnif + demUnif + pcgsp1000 + termLimits + changeForeign, data
## = imm.data)
## weights: nb2listw(state.nb, style = "W")
## 
## LMerr = 1.6814, df = 1, p-value = 0.1947
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = immig0511 ~ pubIdeolCCES + sqrt(squireProfess)
## + repUnif + demUnif + pcgsp1000 + termLimits + changeForeign, data
## = imm.data)
## weights: nb2listw(state.nb, style = "W")
## 
## LMlag = 0.11809, df = 1, p-value = 0.7311
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = immig0511 ~ pubIdeolCCES + sqrt(squireProfess)
## + repUnif + demUnif + pcgsp1000 + termLimits + changeForeign, data
## = imm.data)
## weights: nb2listw(state.nb, style = "W")
## 
## RLMerr = 2.1739, df = 1, p-value = 0.1404
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = immig0511 ~ pubIdeolCCES + sqrt(squireProfess)
## + repUnif + demUnif + pcgsp1000 + termLimits + changeForeign, data
## = imm.data)
## weights: nb2listw(state.nb, style = "W")
## 
## RLMlag = 0.61066, df = 1, p-value = 0.4345
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = immig0511 ~ pubIdeolCCES + sqrt(squireProfess)
## + repUnif + demUnif + pcgsp1000 + termLimits + changeForeign, data
## = imm.data)
## weights: nb2listw(state.nb, style = "W")
## 
## SARMA = 2.292, df = 2, p-value = 0.3179
#Fit Conditional Autoregressive Error Model
hisp.car<-spautolm(immig0511~pubIdeolCCES+sqrt(squireProfess)+repUnif+demUnif+
                     pcgsp1000+termLimits+changeForeign, data=imm.data, 
                   listw=nb2listw(state.nb, style='C'), family='CAR')
summary(hisp.car)
## 
## Call: spautolm(formula = immig0511 ~ pubIdeolCCES + sqrt(squireProfess) + 
##     repUnif + demUnif + pcgsp1000 + termLimits + changeForeign, 
##     data = imm.data, listw = nb2listw(state.nb, style = "C"), 
##     family = "CAR")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.25275 -0.36994  0.10529  0.34431  0.79470 
## 
## Coefficients: 
##                       Estimate Std. Error z value  Pr(>|z|)
## (Intercept)         -1.3881630  0.6414357 -2.1642 0.0304528
## pubIdeolCCES         0.0311108  0.0111077  2.8008 0.0050971
## sqrt(squireProfess)  2.0038527  0.6392588  3.1347 0.0017206
## repUnif              0.0396493  0.0339832  1.1667 0.2433185
## demUnif              0.0313643  0.0362122  0.8661 0.3864211
## pcgsp1000            0.0315273  0.0102941  3.0627 0.0021939
## termLimits          -0.2702251  0.1643909 -1.6438 0.1002184
## changeForeign       -0.0045352  0.0012220 -3.7113 0.0002062
## 
## Lambda: -0.65412 LR test value: 1.4668 p-value: 0.22586 
## Numerical Hessian standard error of lambda: 0.50925 
## 
## Log likelihood: -36.40544 
## ML residual variance (sigma squared): 0.25519, (sigma: 0.50516)
## Number of observations: 48 
## Number of parameters estimated: 10 
## AIC: 92.811
#Fit Simultaneous Autoregressive Error Model
hisp.sar<-spautolm(immig0511~pubIdeolCCES+sqrt(squireProfess)+repUnif+demUnif+
                     pcgsp1000+termLimits+changeForeign, data=imm.data, 
                   listw=nb2listw(state.nb, style='C'), family='SAR')
summary(hisp.sar)
## 
## Call: spautolm(formula = immig0511 ~ pubIdeolCCES + sqrt(squireProfess) + 
##     repUnif + demUnif + pcgsp1000 + termLimits + changeForeign, 
##     data = imm.data, listw = nb2listw(state.nb, style = "C"), 
##     family = "SAR")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.283042 -0.276762  0.097122  0.362259  0.827288 
## 
## Coefficients: 
##                       Estimate Std. Error z value   Pr(>|z|)
## (Intercept)         -1.4337314  0.6344365 -2.2599  0.0238306
## pubIdeolCCES         0.0292941  0.0109771  2.6687  0.0076155
## sqrt(squireProfess)  2.1195678  0.6203145  3.4169  0.0006333
## repUnif              0.0393386  0.0337906  1.1642  0.2443480
## demUnif              0.0342368  0.0361459  0.9472  0.3435454
## pcgsp1000            0.0310867  0.0101159  3.0730  0.0021189
## termLimits          -0.2778785  0.1608923 -1.7271  0.0841482
## changeForeign       -0.0046582  0.0011823 -3.9397 0.00008157
## 
## Lambda: -0.42474 LR test value: 1.8893 p-value: 0.16928 
## Numerical Hessian standard error of lambda: 0.30305 
## 
## Log likelihood: -36.19418 
## ML residual variance (sigma squared): 0.25469, (sigma: 0.50467)
## Number of observations: 48 
## Number of parameters estimated: 10 
## AIC: 92.388
#Fit spatial error model from Ward & Gleditsch
hisp.err<-errorsarlm(immig0511~pubIdeolCCES+sqrt(squireProfess)+repUnif+demUnif+
                       pcgsp1000+termLimits+changeForeign, data=imm.data, 
                     listw=nb2listw(state.nb, style='C'), method='eigen')
summary(hisp.err)
## 
## Call:
## errorsarlm(formula = immig0511 ~ pubIdeolCCES + sqrt(squireProfess) + 
##     repUnif + demUnif + pcgsp1000 + termLimits + changeForeign, 
##     data = imm.data, listw = nb2listw(state.nb, style = "C"), 
##     method = "eigen")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.283042 -0.276762  0.097122  0.362259  0.827288 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                       Estimate Std. Error z value   Pr(>|z|)
## (Intercept)         -1.4337314  0.6344365 -2.2599  0.0238306
## pubIdeolCCES         0.0292941  0.0109771  2.6687  0.0076155
## sqrt(squireProfess)  2.1195679  0.6203145  3.4169  0.0006333
## repUnif              0.0393386  0.0337906  1.1642  0.2443480
## demUnif              0.0342368  0.0361459  0.9472  0.3435454
## pcgsp1000            0.0310867  0.0101159  3.0730  0.0021189
## termLimits          -0.2778785  0.1608923 -1.7271  0.0841482
## changeForeign       -0.0046582  0.0011823 -3.9397 0.00008157
## 
## Lambda: -0.42474, LR test value: 1.8893, p-value: 0.16928
## Asymptotic standard error: 0.22908
##     z-value: -1.8541, p-value: 0.063718
## Wald statistic: 3.4379, p-value: 0.063718
## 
## Log likelihood: -36.19418 for error model
## ML residual variance (sigma squared): 0.25469, (sigma: 0.50467)
## Number of observations: 48 
## Number of parameters estimated: 10 
## AIC: 92.388, (AIC for lm: 92.278)
#Fit spatial lag model from Ward & Gleditsch
hisp.lag<-lagsarlm(immig0511~pubIdeolCCES+sqrt(squireProfess)+repUnif+demUnif+
                     pcgsp1000+termLimits+changeForeign, data=imm.data, 
                   listw=nb2listw(state.nb, style='C'), method='eigen')
summary(hisp.lag)
## 
## Call:lagsarlm(formula = immig0511 ~ pubIdeolCCES + sqrt(squireProfess) + 
##     repUnif + demUnif + pcgsp1000 + termLimits + changeForeign, 
##     data = imm.data, listw = nb2listw(state.nb, style = "C"), 
##     method = "eigen")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.30540 -0.37566  0.06357  0.36865  0.80318 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                       Estimate Std. Error z value  Pr(>|z|)
## (Intercept)         -1.2334590  0.6558271 -1.8808 0.0600034
## pubIdeolCCES         0.0394609  0.0118007  3.3440 0.0008259
## sqrt(squireProfess)  1.8376377  0.7026897  2.6151 0.0089189
## repUnif              0.0384534  0.0338370  1.1364 0.2557769
## demUnif              0.0232128  0.0361121  0.6428 0.5203560
## pcgsp1000            0.0326916  0.0106907  3.0580 0.0022285
## termLimits          -0.2932435  0.1780584 -1.6469 0.0995795
## changeForeign       -0.0045720  0.0013929 -3.2823 0.0010297
## 
## Rho: -0.14321, LR test value: 0.7304, p-value: 0.39275
## Asymptotic standard error: 0.15683
##     z-value: -0.91314, p-value: 0.36117
## Wald statistic: 0.83382, p-value: 0.36117
## 
## Log likelihood: -36.77362 for lag model
## ML residual variance (sigma squared): 0.26979, (sigma: 0.51942)
## Number of observations: 48 
## Number of parameters estimated: 10 
## AIC: 93.547, (AIC for lm: 92.278)
## LM test for residual autocorrelation
## test value: 0.83283, p-value: 0.36146

0.0.6 Hierarchical Spatial Modeling (TODO)

https://aheblog.com/2016/12/07/geostatistical-modelling-with-r-and-stan/

Hierarchical Modeling and Analysis for Spatial Data, 2nd ed. (ISBN-13: 978-1-4398-1917-3), by S. Banerjee, B.P. Carlin and A.E. Gelfand, Boca Raton, FL: Chapman and Hall/CRC Press, 2015.

https://www.counterpointstat.com/hierarchical-modeling-and-analysis-for-spatial-data.html/

0.0.7 Spatial Misalignment

###Replication of Scallop example in R from chapter 2 of Banerjee, Carlin, & Gelfand 2004###

rm(list=ls())

#Create random observations within the block to krige
L <- 100
lat.min <- 39.5
lat.max <- 40.0
long.min <- -73.0
long.max <- -72.5
lat.vec <- runif(L,lat.min,lat.max)
long.vec <- runif(L,long.min,long.max)

#input data
myscallops <- read.table("myscallops.txt", header=T)
#myscallops <- read.table(file.choose(),header=T)

#load libraries
library(akima) #for drawing image/surface plots
library(geoR) #for variogram analyses

#create an image plot with contours #SKIP
int.scp <- interp(myscallops$long, myscallops$lat, myscallops$lgcatch)
image(int.scp, xlim=range(myscallops$long), ylim=range(myscallops$lat))
contour(int.scp, add=T)
polygon(x=c(-72.5,-72.5,-73,-73), y=c(40,39.5,39.5,40),border='cyan')

#What would you guess the average is?

#create a surface plot #SKIP
persp(int.scp, xlim=range(myscallops$long), ylim=range(myscallops$lat))

#define a "geodata" object
obj <- cbind(myscallops$long, myscallops$lat, myscallops$lgcatch)
scallops.geo <- as.geodata(obj, coords.col=1:2, data.col=3)

#create an empirical semivariogram
scallops.var <- variog(scallops.geo, estimator.type="classical")
## variog: computing omnidirectional variogram
scallops.var
## $u
##  [1] 0.1115405 0.3346215 0.5577025 0.7807835 1.0038645 1.2269455 1.4500265
##  [8] 1.6731075 1.8961885 2.1192695 2.3423505 2.5654315 2.7885125
## 
## $v
##  [1] 3.048744 4.361375 4.520836 4.924096 5.697811 5.175957 4.888780
##  [8] 4.345750 4.455269 4.224356 4.460676 5.338326 7.391131
## 
## $n
##  [1]  669 1538 1566 1534 1410 1267 1011  724  465  329  218  109   37
## 
## $sd
##  [1] 5.736655 6.292244 5.600070 6.288254 7.001918 6.329819 6.570810
##  [8] 5.975687 5.815976 4.899078 5.451433 5.928083 7.090025
## 
## $bins.lim
##  [1] 0.000000000001 0.223081002026 0.446162004053 0.669243006079
##  [5] 0.892324008105 1.115405010132 1.338486012158 1.561567014185
##  [9] 1.784648016211 2.007729018237 2.230810020264 2.453891022290
## [13] 2.676972024316 2.900053026343
## 
## $ind.bin
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 
## $var.mark
## [1] 4.720595
## 
## $beta.ols
## [1] 3.482623
## 
## $output.type
## [1] "bin"
## 
## $max.dist
## [1] 2.900053
## 
## $estimator.type
## [1] "classical"
## 
## $n.data
## [1] 148
## 
## $lambda
## [1] 1
## 
## $trend
## [1] "cte"
## 
## $pairs.min
## [1] 2
## 
## $nugget.tolerance
## [1] 1e-12
## 
## $direction
## [1] "omnidirectional"
## 
## $tolerance
## [1] "none"
## 
## $uvec
##  [1] 0.1115405 0.3346215 0.5577025 0.7807835 1.0038645 1.2269455 1.4500265
##  [8] 1.6731075 1.8961885 2.1192695 2.3423505 2.5654315 2.7885125
## 
## $call
## variog(geodata = scallops.geo, estimator.type = "classical")
## 
## attr(,"class")
## [1] "variogram"
#create a robust empirical semivariogram
scallops.var.robust <- variog(scallops.geo, estimator.type="modulus")
## variog: computing omnidirectional variogram
#plot the two forms of the semivariogram
par(mfrow=c(2,1)) #align separate plots on the same space, two rows by one column
plot(scallops.var)
plot(scallops.var.robust)

#Fit a parametric variogram model (assume an exponential model) using weighted least squares
scallops.var.fit <- variofit(scallops.var.robust, ini.cov.pars = c(1.0,2.0),
                             cov.model="exponential", fix.nugget=FALSE, nugget=1.0)
## variofit: covariance model used is exponential 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim
## Warning in variofit(scallops.var.robust, ini.cov.pars = c(1, 2), cov.model
## = "exponential", : unreasonable initial value for sigmasq + nugget (too
## low)
#Fit the exponential model using MLE
scallops.lik.fit <- likfit(scallops.geo, ini.cov.pars=c(1.0,2.0), 
                           cov.model = "exponential", trend = "cte", 
                           fix.nugget = FALSE, nugget = 1.0, 
                           nospatial = TRUE, lik.method = "ML") 
## kappa not used for the exponential correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.
#note: "lik.method" is a correction of a typo from the text
scallops.lik.fit
## likfit: estimated model parameters:
##     beta    tausq  sigmasq      phi 
## "2.3748" "0.0947" "5.7675" "0.2338" 
## Practical Range with cor=0.05 for asymptotic range: 0.7005547
## 
## likfit: maximised log-likelihood = -285.9
#CLASSICAL KRIGING
#First we must create locations to predict#
#Source: http://www.biostat.umn.edu/~sudiptob/Software/Book.R
u.loci <- cbind(long.vec, lat.vec)
scallops.classical <- krige.conv(geodata = scallops.geo, 
                                 locations = u.loci, 
                                 borders = NULL, 
                                 krige = krige.control(cov.pars=c(5.7675,0.2338))) #cov.pars=c(sigmasq,phi)
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
scallops.classical$predict
##   [1] 5.6419341 1.7432806 6.9119948 2.7666652 6.5711829 4.1043280 0.7977967
##   [8] 5.8868934 4.0250767 5.7365395 6.6456492 6.8274057 7.1243194 4.4532750
##  [15] 6.1682690 6.9426476 3.9297373 3.8282980 6.1188330 6.3802473 6.7251522
##  [22] 2.4312663 7.3685875 6.3593291 0.5356865 4.9771101 6.6788021 0.7907546
##  [29] 5.6073281 5.9049236 4.1623659 4.4916631 6.3101120 5.4695400 6.6981826
##  [36] 7.2099658 4.7542212 2.4198100 6.5279612 1.2411826 4.6001424 2.4825210
##  [43] 5.6267846 3.5800146 3.1160066 4.0481120 3.1276008 7.6949181 6.5317052
##  [50] 3.4223027 4.5079967 7.0502903 6.6481609 6.5888756 7.3429442 4.3821784
##  [57] 2.3024721 4.8568017 5.2744735 5.5499560 1.8034411 5.2290168 7.1828775
##  [64] 7.8271450 0.8225617 7.6170001 6.4242428 7.3011851 0.5264109 2.8854938
##  [71] 2.9805537 4.3563635 6.8930906 7.2032837 6.3981411 7.0144640 2.4763012
##  [78] 3.8205222 4.0216234 3.4520632 5.1477982 7.2604121 1.0833592 3.6994653
##  [85] 3.6938089 5.3850278 4.4454305 4.7055334 5.8076492 3.6998101 7.0274229
##  [92] 5.1036242 4.9551615 5.9113040 0.7688712 1.8666832 0.8744410 4.7890868
##  [99] 3.7841463 7.4429925
mean(scallops.classical$predict) #expected average scallop catch in the block
## [1] 4.816904
quantile(scallops.classical$predict, c(.5,.05,.95)) #median and 90% credible interval
##       50%        5%       95% 
## 5.0403671 0.8213234 7.3442264
#Bayesian Kriging Model--Estimation and Prediction
scallops.bayes2 <- krige.bayes(scallops.geo,
    locations = u.loci, #defaults to locations="no" 
    borders = NULL,
    model = model.control(trend.d = "cte",
        cov.model = "exponential"),
    prior = prior.control(beta.prior = "flat",
        sigmasq.prior = "reciprocal",
        tausq.rel.prior = "uniform",
        tausq.rel.discrete = seq(from=0.0, to=1.0, by=0.01)))
## krige.bayes: model with constant mean
## krige.bayes: computing the discrete posterior of phi/tausq.rel
## krige.bayes: argument `phi.discrete` not provided, using default values
## krige.bayes: computing the posterior probabilities.
##              Number of parameter sets:  5050 
## 1, 101, 201, 301, 401, 501, 601, 701, 801, 901, 1001, 1101, 1201, 1301, 1401, 1501, 1601, 1701, 1801, 1901, 2001, 2101, 2201, 2301, 2401, 2501, 2601, 2701, 2801, 2901, 3001, 3101, 3201, 3301, 3401, 3501, 3601, 3701, 3801, 3901, 4001, 4101, 4201, 4301, 4401, 4501, 4601, 4701, 4801, 4901, 5001, 
## 
## krige.bayes: sampling from posterior distribution
## krige.bayes: sample from the (joint) posterior of phi and tausq.rel
##                [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
## phi       0.1160021 0.2320042 0.3480064 0.4640085 0.5800106 0.6960127
## tausq.rel 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## frequency 1.0000000 9.0000000 9.0000000 9.0000000 6.0000000 6.0000000
##                [,7]     [,8]     [,9]    [,10]    [,11]    [,12]    [,13]
## phi       0.8120148 0.928017 1.044019 1.160021 1.276023 1.392025 1.508028
## tausq.rel 0.0000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
## frequency 8.0000000 7.000000 2.000000 3.000000 4.000000 2.000000 1.000000
##             [,14]    [,15]    [,16]    [,17]    [,18]   [,19]    [,20]
## phi       1.62403 1.740032 1.856034 1.972036 2.088038 2.20404 2.668049
## tausq.rel 0.00000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000
## frequency 4.00000 3.000000 1.000000 7.000000 2.000000 1.00000 2.000000
##              [,21]    [,22]    [,23]    [,24]    [,25]    [,26]    [,27]
## phi       2.784051 3.016055 3.132057 3.248059 3.364062 3.480064 3.712068
## tausq.rel 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
## frequency 4.000000 6.000000 2.000000 4.000000 2.000000 1.000000 3.000000
##             [,28]    [,29]    [,30]    [,31]    [,32]    [,33]    [,34]
## phi       3.82807 3.944072 4.060074 4.176076 4.408081 4.524083 4.640085
## tausq.rel 0.00000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
## frequency 1.00000 1.000000 2.000000 3.000000 1.000000 5.000000 2.000000
##              [,35]    [,36]    [,37]    [,38]  [,39]    [,40]      [,41]
## phi       4.756087 4.988091 5.104093 5.336098 5.4521 5.568102  0.2320042
## tausq.rel 0.000000 0.000000 0.000000 0.000000 0.0000 0.000000  0.0100000
## frequency 3.000000 1.000000 1.000000 2.000000 2.0000 2.000000 17.0000000
##                [,42]      [,43]     [,44]     [,45]      [,46]     [,47]
## phi        0.3480064  0.4640085 0.5800106 0.6960127  0.8120148  0.928017
## tausq.rel  0.0100000  0.0100000 0.0100000 0.0100000  0.0100000  0.010000
## frequency 17.0000000 10.0000000 7.0000000 4.0000000 11.0000000 11.000000
##              [,48]    [,49]    [,50]    [,51]    [,52]   [,53]    [,54]
## phi       1.044019 1.160021 1.276023 1.392025 1.508028 1.62403 1.740032
## tausq.rel 0.010000 0.010000 0.010000 0.010000 0.010000 0.01000 0.010000
## frequency 5.000000 7.000000 6.000000 3.000000 4.000000 4.00000 4.000000
##              [,55]    [,56]    [,57]   [,58]    [,59]    [,60]    [,61]
## phi       1.856034 1.972036 2.088038 2.20404 2.320042 2.436045 2.552047
## tausq.rel 0.010000 0.010000 0.010000 0.01000 0.010000 0.010000 0.010000
## frequency 4.000000 4.000000 4.000000 7.00000 7.000000 5.000000 3.000000
##              [,62]    [,63]    [,64]    [,65]    [,66]    [,67]    [,68]
## phi       2.668049 2.784051 2.900053 3.016055 3.132057 3.248059 3.364062
## tausq.rel 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000
## frequency 5.000000 3.000000 2.000000 5.000000 1.000000 2.000000 2.000000
##              [,69]    [,70]   [,71]    [,72]    [,73]    [,74]    [,75]
## phi       3.480064 3.596066 3.82807 3.944072 4.060074 4.176076 4.292078
## tausq.rel 0.010000 0.010000 0.01000 0.010000 0.010000 0.010000 0.010000
## frequency 5.000000 2.000000 1.00000 3.000000 5.000000 5.000000 2.000000
##              [,76]    [,77]    [,78]    [,79]    [,80]    [,81]    [,82]
## phi       4.408081 4.640085 4.756087 4.872089 4.988091 5.104093 5.336098
## tausq.rel 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000
## frequency 2.000000 1.000000 2.000000 3.000000 4.000000 2.000000 1.000000
##              [,83]    [,84]      [,85]      [,86]      [,87]      [,88]
## phi       5.568102 5.684104  0.2320042  0.3480064  0.4640085  0.5800106
## tausq.rel 0.010000 0.010000  0.0200000  0.0200000  0.0200000  0.0200000
## frequency 1.000000 2.000000 17.0000000 10.0000000 14.0000000 15.0000000
##                [,89]      [,90]    [,91]    [,92]    [,93]    [,94]
## phi        0.6960127  0.8120148 0.928017 1.044019 1.160021 1.276023
## tausq.rel  0.0200000  0.0200000 0.020000 0.020000 0.020000 0.020000
## frequency 11.0000000 10.0000000 4.000000 3.000000 6.000000 7.000000
##              [,95]    [,96]   [,97]    [,98]    [,99]   [,100]   [,101]
## phi       1.392025 1.508028 1.62403 1.740032 1.856034 1.972036 2.088038
## tausq.rel 0.020000 0.020000 0.02000 0.020000 0.020000 0.020000 0.020000
## frequency 2.000000 6.000000 4.00000 2.000000 5.000000 4.000000 5.000000
##            [,102]   [,103]   [,104]   [,105]   [,106]   [,107]   [,108]
## phi       2.20404 2.320042 2.436045 2.552047 2.668049 2.784051 2.900053
## tausq.rel 0.02000 0.020000 0.020000 0.020000 0.020000 0.020000 0.020000
## frequency 2.00000 1.000000 4.000000 4.000000 2.000000 1.000000 1.000000
##             [,109]   [,110]   [,111]   [,112]   [,113]   [,114]   [,115]
## phi       3.016055 3.248059 3.596066 3.712068 4.060074 4.408081 4.524083
## tausq.rel 0.020000 0.020000 0.020000 0.020000 0.020000 0.020000 0.020000
## frequency 1.000000 1.000000 1.000000 3.000000 1.000000 1.000000 1.000000
##             [,116]     [,117]     [,118]     [,119]    [,120]    [,121]
## phi       4.640085  0.2320042  0.3480064  0.4640085 0.5800106 0.6960127
## tausq.rel 0.020000  0.0300000  0.0300000  0.0300000 0.0300000 0.0300000
## frequency 1.000000 14.0000000 18.0000000 13.0000000 7.0000000 6.0000000
##               [,122]   [,123]   [,124]   [,125]   [,126]   [,127]   [,128]
## phi        0.8120148 0.928017 1.044019 1.160021 1.276023 1.392025 1.508028
## tausq.rel  0.0300000 0.030000 0.030000 0.030000 0.030000 0.030000 0.030000
## frequency 10.0000000 6.000000 4.000000 5.000000 3.000000 2.000000 4.000000
##            [,129]   [,130]   [,131]   [,132]  [,133]   [,134]     [,135]
## phi       1.62403 1.740032 1.856034 1.972036 2.20404 4.756087  0.2320042
## tausq.rel 0.03000 0.030000 0.030000 0.030000 0.03000 0.030000  0.0400000
## frequency 1.00000 2.000000 1.000000 4.000000 1.00000 1.000000 11.0000000
##               [,136]    [,137]    [,138]    [,139]    [,140]   [,141]
## phi        0.3480064 0.4640085 0.5800106 0.6960127 0.8120148 0.928017
## tausq.rel  0.0400000 0.0400000 0.0400000 0.0400000 0.0400000 0.040000
## frequency 11.0000000 6.0000000 9.0000000 7.0000000 6.0000000 9.000000
##             [,142]   [,143]   [,144]   [,145]   [,146]  [,147]   [,148]
## phi       1.044019 1.160021 1.276023 1.392025 1.508028 1.62403 1.740032
## tausq.rel 0.040000 0.040000 0.040000 0.040000 0.040000 0.04000 0.040000
## frequency 5.000000 3.000000 5.000000 1.000000 4.000000 1.00000 1.000000
##             [,149]   [,150]   [,151]   [,152]     [,153]    [,154]
## phi       2.436045 3.016055 3.248059 3.712068  0.2320042 0.3480064
## tausq.rel 0.040000 0.040000 0.040000 0.040000  0.0500000 0.0500000
## frequency 1.000000 1.000000 1.000000 1.000000 10.0000000 8.0000000
##              [,155]     [,156]    [,157]    [,158]   [,159]   [,160]
## phi       0.4640085  0.5800106 0.6960127 0.8120148 0.928017 1.044019
## tausq.rel 0.0500000  0.0500000 0.0500000 0.0500000 0.050000 0.050000
## frequency 9.0000000 11.0000000 6.0000000 3.0000000 3.000000 2.000000
##             [,161]   [,162]   [,163]   [,164]    [,165]     [,166]
## phi       1.160021 1.276023 1.392025 1.740032 0.2320042  0.3480064
## tausq.rel 0.050000 0.050000 0.050000 0.050000 0.0600000  0.0600000
## frequency 1.000000 3.000000 2.000000 1.000000 8.0000000 13.0000000
##              [,167]    [,168]    [,169]    [,170]   [,171]   [,172]
## phi       0.4640085 0.5800106 0.6960127 0.8120148 0.928017 1.044019
## tausq.rel 0.0600000 0.0600000 0.0600000 0.0600000 0.060000 0.060000
## frequency 7.0000000 7.0000000 6.0000000 5.0000000 3.000000 1.000000
##             [,173]   [,174]   [,175]   [,176]     [,177]     [,178]
## phi       1.160021 1.276023 1.392025 1.856034  0.2320042  0.3480064
## tausq.rel 0.060000 0.060000 0.060000 0.060000  0.0700000  0.0700000
## frequency 1.000000 1.000000 2.000000 1.000000 11.0000000 12.0000000
##              [,179]    [,180]    [,181]    [,182]   [,183]   [,184]
## phi       0.4640085 0.5800106 0.6960127 0.8120148 0.928017 1.044019
## tausq.rel 0.0700000 0.0700000 0.0700000 0.0700000 0.070000 0.070000
## frequency 7.0000000 5.0000000 5.0000000 1.0000000 1.000000 3.000000
##             [,185]   [,186]    [,187]    [,188]    [,189]    [,190]
## phi       1.392025 1.508028 0.1160021 0.2320042 0.3480064 0.4640085
## tausq.rel 0.070000 0.070000 0.0800000 0.0800000 0.0800000 0.0800000
## frequency 1.000000 1.000000 1.0000000 5.0000000 8.0000000 8.0000000
##              [,191]    [,192]    [,193]   [,194]    [,195]    [,196]
## phi       0.5800106 0.6960127 0.8120148 1.276023 0.2320042 0.3480064
## tausq.rel 0.0800000 0.0800000 0.0800000 0.080000 0.0900000 0.0900000
## frequency 3.0000000 1.0000000 2.0000000 2.000000 1.0000000 9.0000000
##              [,197]    [,198]    [,199]   [,200]   [,201]    [,202]
## phi       0.4640085 0.5800106 0.8120148 0.928017 1.740032 0.2320042
## tausq.rel 0.0900000 0.0900000 0.0900000 0.090000 0.090000 0.1000000
## frequency 2.0000000 4.0000000 1.0000000 2.000000 1.000000 7.0000000
##              [,203]    [,204]    [,205]    [,206]    [,207]    [,208]
## phi       0.3480064 0.4640085 0.5800106 0.8120148 0.2320042 0.3480064
## tausq.rel 0.1000000 0.1000000 0.1000000 0.1000000 0.1100000 0.1100000
## frequency 4.0000000 1.0000000 2.0000000 1.0000000 4.0000000 5.0000000
##              [,209]    [,210]    [,211]    [,212]   [,213]    [,214]
## phi       0.4640085 0.5800106 0.6960127 0.8120148 0.928017 0.2320042
## tausq.rel 0.1100000 0.1100000 0.1100000 0.1100000 0.110000 0.1200000
## frequency 5.0000000 1.0000000 1.0000000 1.0000000 1.000000 4.0000000
##              [,215]    [,216]    [,217]    [,218]   [,219]   [,220]
## phi       0.3480064 0.4640085 0.5800106 0.8120148 0.928017 1.044019
## tausq.rel 0.1200000 0.1200000 0.1200000 0.1200000 0.120000 0.120000
## frequency 3.0000000 3.0000000 1.0000000 1.0000000 1.000000 1.000000
##             [,221]    [,222]    [,223]    [,224]    [,225]    [,226]
## phi       1.160021 0.2320042 0.3480064 0.4640085 0.5800106 0.6960127
## tausq.rel 0.120000 0.1300000 0.1300000 0.1300000 0.1300000 0.1300000
## frequency 1.000000 3.0000000 3.0000000 2.0000000 1.0000000 2.0000000
##              [,227]    [,228]    [,229]    [,230]    [,231]    [,232]
## phi       0.8120148 0.2320042 0.3480064 0.4640085 0.5800106 0.6960127
## tausq.rel 0.1300000 0.1400000 0.1400000 0.1400000 0.1400000 0.1400000
## frequency 1.0000000 2.0000000 3.0000000 4.0000000 1.0000000 1.0000000
##              [,233]   [,234]    [,235]    [,236]    [,237]    [,238]
## phi       0.8120148 0.928017 0.2320042 0.3480064 0.4640085 0.5800106
## tausq.rel 0.1400000 0.140000 0.1500000 0.1500000 0.1500000 0.1500000
## frequency 1.0000000 1.000000 5.0000000 1.0000000 3.0000000 1.0000000
##              [,239]    [,240]    [,241]    [,242]    [,243]    [,244]
## phi       0.6960127 0.2320042 0.3480064 0.4640085 0.2320042 0.3480064
## tausq.rel 0.1500000 0.1600000 0.1600000 0.1600000 0.1700000 0.1700000
## frequency 1.0000000 2.0000000 1.0000000 2.0000000 4.0000000 1.0000000
##              [,245]    [,246]    [,247]    [,248]    [,249]    [,250]
## phi       0.4640085 0.2320042 0.3480064 0.2320042 0.3480064 0.2320042
## tausq.rel 0.1700000 0.1800000 0.1800000 0.1900000 0.1900000 0.2000000
## frequency 1.0000000 2.0000000 2.0000000 2.0000000 1.0000000 1.0000000
##              [,251]    [,252]    [,253]    [,254]    [,255]    [,256]
## phi       0.3480064 0.2320042 0.3480064 0.2320042 0.3480064 0.3480064
## tausq.rel 0.2000000 0.2200000 0.2200000 0.2300000 0.2300000 0.2400000
## frequency 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
##              [,257]    [,258]    [,259]    [,260]    [,261]
## phi       0.2320042 0.3480064 0.2320042 0.3480064 0.2320042
## tausq.rel 0.2800000 0.2800000 0.2900000 0.3200000 0.5400000
## frequency 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## 
## krige.bayes: starting prediction at the provided locations
## krige.bayes: phi/tausq.rel samples for the predictive are same as for the posterior 
## krige.bayes: computing moments of the predictive distribution
## krige.bayes: sampling from the predictive
##              Number of parameter sets:  261 
## 1, 11, 21, 31, 41, 51, 61, 71, 81, 91, 101, 111, 121, 131, 141, 151, 161, 171, 181, 191, 201, 211, 221, 231, 241, 251, 261, 
## krige.bayes: preparing summaries of the predictive distribution
#TURN TO RESULTS PAGE

# Projected number of catches from Bayesian model
mean(scallops.bayes2$predictive$mean)
## [1] 4.811959
quantile(scallops.bayes2$predictive$mean, c(.5,.05,.95))
##      50%       5%      95% 
## 5.024128 1.248310 7.123376
#Form the quantiles to interpret the model results
out <- scallops.bayes2$posterior
out <- out$sample
beta.qnt <- quantile(out$beta, c(0.50,0.025,0.975))
phi.qnt <-quantile(out$phi, c(0.50,0.025,0.975))
sigmasq.qnt <- quantile(out$sigmasq, c(0.50,0.025,0.975))
tausq.rel.qnt <- quantile(out$tausq.rel, c(0.50, 0.025, 0.975))
beta.qnt
##       50%      2.5%     97.5% 
##  1.766796 -7.892588  7.216852
phi.qnt
##       50%      2.5%     97.5% 
## 0.6960127 0.2320042 4.7560870
sigmasq.qnt
##       50%      2.5%     97.5% 
## 11.966236  4.258111 86.864834
tausq.rel.qnt
##     50%    2.5%   97.5% 
## 0.03000 0.00000 0.16025

0.0.8 Spatiotemporal Modeling

###Replication of:
#Franzese and Hays. 2007. "Spatial Econometric Models of Cross-Sectional 
#Interdependence in Political Science Panel and Time-Series-Cross-Section Data."
#Political Analysis 15:140-164.
###Usage of these data or code MUST cite the original paper.


#clean up
rm(list=ls())

#preferences
options(scipen=8)

#load required libraries
library(spdep) 
library(maps)
library(foreign)
library(rgdal)
library(plm)
## Loading required package: Formula
#Load Data
fh<-read.csv('almMLire.csv',header=TRUE)

#Load Adjacency Matrix
adj<-as.matrix(read.csv('almWeights.csv',header=FALSE))

#Standardize Adjacency Matrix
noNeigh<-apply(adj,1,sum)
noNeigh[7]<-.01 #no neighbors for Greece, can't have 0 in the denominator
wmat <- adj/noNeigh

#Spatial Lag of the Outcome Variable, Using Kronecker Product
N<-15
T<-12
ident<-diag(1,nrow=T)
bigW<-(ident%x%wmat)
fh$SpatLag<-(ident%x%wmat)%*%fh$lnlmtue

#####REPLICATE TABLE 4#####
###First Model: Temporal Lag Only###
#dummied-out
mod.4.1.a<-lm(lnlmtue~lnlmtue_1+DENSITY+DEIND+lngdp_pc+UR+TRADE+FDI+LLVOTE+
                LEFTC+TCDEMC+GOVCON+OLDAGE+as.factor(cc)+as.factor(year),data=fh)
summary(mod.4.1.a)
## 
## Call:
## lm(formula = lnlmtue ~ lnlmtue_1 + DENSITY + DEIND + lngdp_pc + 
##     UR + TRADE + FDI + LLVOTE + LEFTC + TCDEMC + GOVCON + OLDAGE + 
##     as.factor(cc) + as.factor(year), data = fh)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.19165 -0.14680  0.01142  0.12235  1.10187 
## 
## Coefficients:
##                       Estimate Std. Error t value         Pr(>|t|)    
## (Intercept)         -1.3331830  3.5995002  -0.370          0.71165    
## lnlmtue_1            0.5812769  0.0779578   7.456 0.00000000000804 ***
## DENSITY              0.0180919  0.0105857   1.709          0.08962 .  
## DEIND                0.0619381  0.0338505   1.830          0.06938 .  
## lngdp_pc            -0.3283057  0.9457078  -0.347          0.72899    
## UR                  -0.0892693  0.0334585  -2.668          0.00852 ** 
## TRADE                0.0007456  0.0085212   0.088          0.93039    
## FDI                 -0.0228762  0.0249063  -0.918          0.35992    
## LLVOTE              -0.0130444  0.0153701  -0.849          0.39748    
## LEFTC                0.0003791  0.0009957   0.381          0.70395    
## TCDEMC              -0.0025288  0.0035060  -0.721          0.47193    
## GOVCON               0.0359730  0.0379986   0.947          0.34540    
## OLDAGE              -0.0106190  0.0585543  -0.181          0.85635    
## as.factor(cc)2      -0.3317578  0.6096939  -0.544          0.58720    
## as.factor(cc)3      -0.3808918  0.4829377  -0.789          0.43160    
## as.factor(cc)4      -0.3877157  0.4384354  -0.884          0.37802    
## as.factor(cc)5       0.5039088  0.4699241   1.072          0.28540    
## as.factor(cc)6       0.7928225  0.2462236   3.220          0.00159 ** 
## as.factor(cc)7       0.1212794  0.4897753   0.248          0.80478    
## as.factor(cc)8       0.0129930  0.5178942   0.025          0.98002    
## as.factor(cc)9       0.0145534  0.5347751   0.027          0.97833    
## as.factor(cc)10     -0.3887082  0.4065054  -0.956          0.34059    
## as.factor(cc)11     -0.0886912  0.5295007  -0.167          0.86722    
## as.factor(cc)12      0.3670486  0.5824373   0.630          0.52958    
## as.factor(cc)13     -0.6791023  0.5635019  -1.205          0.23015    
## as.factor(cc)14      0.4838704  0.4304817   1.124          0.26290    
## as.factor(cc)15     -0.3977154  0.3402356  -1.169          0.24439    
## as.factor(year)1988  0.1158841  0.1193023   0.971          0.33303    
## as.factor(year)1989  0.3078653  0.1285753   2.394          0.01795 *  
## as.factor(year)1990  0.1756465  0.1395881   1.258          0.21034    
## as.factor(year)1991 -0.0388208  0.1502657  -0.258          0.79651    
## as.factor(year)1992  0.0430456  0.1654094   0.260          0.79506    
## as.factor(year)1993  0.0050429  0.1860829   0.027          0.97842    
## as.factor(year)1994  0.0927372  0.2010415   0.461          0.64530    
## as.factor(year)1995  0.1430783  0.2110666   0.678          0.49895    
## as.factor(year)1996  0.1294835  0.2295124   0.564          0.57353    
## as.factor(year)1997  0.1665405  0.2500523   0.666          0.50648    
## as.factor(year)1998  0.2874073  0.2683866   1.071          0.28604    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3127 on 142 degrees of freedom
## Multiple R-squared:  0.8986, Adjusted R-squared:  0.8722 
## F-statistic: 34.02 on 37 and 142 DF,  p-value: < 2.2e-16
#within estimator
mod.4.1.b<-plm(lnlmtue~lnlmtue_1+DENSITY+DEIND+lngdp_pc+UR+TRADE+FDI+LLVOTE+
                 LEFTC+TCDEMC+GOVCON+OLDAGE+as.factor(year),data=fh,
               index=c("cc","year"),model="within")
summary(mod.4.1.b)
## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = lnlmtue ~ lnlmtue_1 + DENSITY + DEIND + lngdp_pc + 
##     UR + TRADE + FDI + LLVOTE + LEFTC + TCDEMC + GOVCON + OLDAGE + 
##     as.factor(year), data = fh, model = "within", index = c("cc", 
##     "year"))
## 
## Balanced Panel: n = 15, T = 12, N = 180
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -1.191652 -0.146805  0.011423  0.122348  1.101871 
## 
## Coefficients:
##                        Estimate  Std. Error t-value          Pr(>|t|)    
## lnlmtue_1            0.58127691  0.07795781  7.4563 0.000000000008038 ***
## DENSITY              0.01809190  0.01058574  1.7091          0.089620 .  
## DEIND                0.06193812  0.03385047  1.8298          0.069384 .  
## lngdp_pc            -0.32830568  0.94570780 -0.3472          0.728990    
## UR                  -0.08926930  0.03345849 -2.6681          0.008516 ** 
## TRADE                0.00074564  0.00852116  0.0875          0.930393    
## FDI                 -0.02287622  0.02490634 -0.9185          0.359920    
## LLVOTE              -0.01304438  0.01537006 -0.8487          0.397484    
## LEFTC                0.00037912  0.00099568  0.3808          0.703948    
## TCDEMC              -0.00252877  0.00350601 -0.7213          0.471932    
## GOVCON               0.03597300  0.03799859  0.9467          0.345403    
## OLDAGE              -0.01061897  0.05855432 -0.1814          0.856350    
## as.factor(year)1988  0.11588407  0.11930226  0.9713          0.333026    
## as.factor(year)1989  0.30786531  0.12857529  2.3944          0.017951 *  
## as.factor(year)1990  0.17564648  0.13958809  1.2583          0.210341    
## as.factor(year)1991 -0.03882076  0.15026568 -0.2583          0.796513    
## as.factor(year)1992  0.04304559  0.16540942  0.2602          0.795058    
## as.factor(year)1993  0.00504291  0.18608285  0.0271          0.978418    
## as.factor(year)1994  0.09273718  0.20104151  0.4613          0.645301    
## as.factor(year)1995  0.14307828  0.21106656  0.6779          0.498949    
## as.factor(year)1996  0.12948348  0.22951236  0.5642          0.573530    
## as.factor(year)1997  0.16654050  0.25005231  0.6660          0.506477    
## as.factor(year)1998  0.28740728  0.26838665  1.0709          0.286045    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    33.616
## Residual Sum of Squares: 13.881
## R-Squared:      0.58708
## Adj. R-Squared: 0.47949
## F-statistic: 8.77806 on 23 and 142 DF, p-value: < 2.22e-16
###Second Model: OLS with a Spatial Lag###
mod.4.2<-plm(lnlmtue~lnlmtue_1+SpatLag+DENSITY+DEIND+lngdp_pc+UR+TRADE+FDI+
               LLVOTE+LEFTC+TCDEMC+GOVCON+OLDAGE+as.factor(year),data=fh,
             index=c("cc","year"),model="within")
summary(mod.4.2)
## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = lnlmtue ~ lnlmtue_1 + SpatLag + DENSITY + DEIND + 
##     lngdp_pc + UR + TRADE + FDI + LLVOTE + LEFTC + TCDEMC + GOVCON + 
##     OLDAGE + as.factor(year), data = fh, model = "within", index = c("cc", 
##     "year"))
## 
## Balanced Panel: n = 15, T = 12, N = 180
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -1.1573255 -0.1512225  0.0015748  0.1438301  1.0948248 
## 
## Coefficients:
##                          Estimate    Std. Error t-value     Pr(>|t|)    
## lnlmtue_1            0.4354879852  0.0808440326  5.3868 0.0000002931 ***
## SpatLag             -0.4376800180  0.1010789502 -4.3301 0.0000281051 ***
## DENSITY              0.0076330163  0.0102684684  0.7433     0.458510    
## DEIND                0.0402701526  0.0323044963  1.2466     0.214618    
## lngdp_pc            -1.3073185100  0.9198428304 -1.4212     0.157455    
## UR                  -0.0974263443  0.0316012054 -3.0830     0.002466 ** 
## TRADE                0.0000062976  0.0080356504  0.0008     0.999376    
## FDI                 -0.0335940322  0.0236120494 -1.4227     0.157017    
## LLVOTE               0.0043235810  0.0150359161  0.2876     0.774113    
## LEFTC               -0.0001629785  0.0009470531 -0.1721     0.863613    
## TCDEMC              -0.0028703711  0.0033064461 -0.8681     0.386807    
## GOVCON               0.0224772985  0.0359607773  0.6251     0.532948    
## OLDAGE               0.0516600127  0.0570484560  0.9055     0.366721    
## as.factor(year)1988  0.2148278697  0.1147769654  1.8717     0.063320 .  
## as.factor(year)1989  0.5480530853  0.1333104306  4.1111 0.0000665474 ***
## as.factor(year)1990  0.5016785281  0.1516218532  3.3087     0.001190 ** 
## as.factor(year)1991  0.2465674015  0.1562525590  1.5780     0.116805    
## as.factor(year)1992  0.3011017340  0.1669491208  1.8036     0.073436 .  
## as.factor(year)1993  0.2106993953  0.1817559899  1.1592     0.248316    
## as.factor(year)1994  0.3030077992  0.1956656343  1.5486     0.123720    
## as.factor(year)1995  0.3742281062  0.2060314733  1.8164     0.071439 .  
## as.factor(year)1996  0.3771705533  0.2238195069  1.6852     0.094171 .  
## as.factor(year)1997  0.4631084695  0.2454991847  1.8864     0.061298 .  
## as.factor(year)1998  0.6604278507  0.2672999375  2.4707     0.014676 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    33.616
## Residual Sum of Squares: 12.252
## R-Squared:      0.63555
## Adj. R-Squared: 0.53733
## F-statistic: 10.2451 on 24 and 141 DF, p-value: < 2.22e-16
###Third Model: Instrumenting the Spatial Lag in Two-Stage Least Squares###
#Spatially Lag the Xs
fh$l.DENSITY<-(ident%x%wmat)%*%fh$DENSITY
fh$l.DEIND<-(ident%x%wmat)%*%fh$DEIND
fh$l.lngdp_pc<-(ident%x%wmat)%*%fh$lngdp_pc
fh$l.UR<-(ident%x%wmat)%*%fh$UR
fh$l.TRADE<-(ident%x%wmat)%*%fh$TRADE
fh$l.FDI<-(ident%x%wmat)%*%fh$FDI
fh$l.LLVOTE<-(ident%x%wmat)%*%fh$LLVOTE
fh$l.LEFTC<-(ident%x%wmat)%*%fh$LEFTC
fh$l.TCDEMC<-(ident%x%wmat)%*%fh$TCDEMC
fh$l.GOVCON<-(ident%x%wmat)%*%fh$GOVCON
fh$l.OLDAGE<-(ident%x%wmat)%*%fh$OLDAGE
fh$l.lnlmtue_1<-(ident%x%wmat)%*%fh$lnlmtue_1

#Model the spatial lag of Y
#2SLS actually calls for modeling all the Xs! Hence, I'm getting somewhat different results.
stage.1<-lm(SpatLag~l.lnlmtue_1+l.DENSITY+l.DEIND+l.lngdp_pc+l.UR+l.TRADE+
              l.FDI+l.LLVOTE+l.LEFTC+l.TCDEMC+l.GOVCON+l.OLDAGE, data=fh)#; summary(stage.1)#######
fh$iv.Lag<-stage.1$fitted.values

#Second Stage Equation
mod.4.3<-plm(lnlmtue~lnlmtue_1+iv.Lag+DENSITY+DEIND+lngdp_pc+UR+TRADE+FDI+
               LLVOTE+LEFTC+TCDEMC+GOVCON+OLDAGE+as.factor(year),data=fh,
             index=c("cc","year"),model="within")
summary(mod.4.3)
## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = lnlmtue ~ lnlmtue_1 + iv.Lag + DENSITY + DEIND + 
##     lngdp_pc + UR + TRADE + FDI + LLVOTE + LEFTC + TCDEMC + GOVCON + 
##     OLDAGE + as.factor(year), data = fh, model = "within", index = c("cc", 
##     "year"))
## 
## Balanced Panel: n = 15, T = 12, N = 180
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -1.2513435 -0.1325319  0.0022381  0.1336672  1.1433427 
## 
## Coefficients:
##                        Estimate  Std. Error t-value      Pr(>|t|)    
## lnlmtue_1            0.49183370  0.08552898  5.7505 0.00000005301 ***
## iv.Lag              -0.29385894  0.12417415 -2.3665      0.019317 *  
## DENSITY              0.01133962  0.01080197  1.0498      0.295619    
## DEIND                0.04422028  0.03414603  1.2950      0.197425    
## lngdp_pc            -0.74604335  0.94734235 -0.7875      0.432304    
## UR                  -0.09255700  0.03295863 -2.8083      0.005688 ** 
## TRADE                0.00295506  0.00843820  0.3502      0.726711    
## FDI                 -0.02923424  0.02465924 -1.1855      0.237803    
## LLVOTE              -0.00222941  0.01580224 -0.1411      0.888007    
## LEFTC                0.00017133  0.00098386  0.1741      0.862006    
## TCDEMC              -0.00294343  0.00345501 -0.8519      0.395698    
## GOVCON               0.03049498  0.03746921  0.8139      0.417093    
## OLDAGE               0.02949127  0.06006905  0.4910      0.624221    
## as.factor(year)1988  0.16224371  0.11903848  1.3630      0.175070    
## as.factor(year)1989  0.40823609  0.13346051  3.0589      0.002660 ** 
## as.factor(year)1990  0.37013026  0.16008512  2.3121      0.022221 *  
## as.factor(year)1991  0.19799317  0.17856377  1.1088      0.269400    
## as.factor(year)1992  0.22534055  0.18009864  1.2512      0.212931    
## as.factor(year)1993  0.17292861  0.19640025  0.8805      0.380092    
## as.factor(year)1994  0.21630155  0.20463542  1.0570      0.292315    
## as.factor(year)1995  0.28224374  0.21589190  1.3073      0.193226    
## as.factor(year)1996  0.26475533  0.23300283  1.1363      0.257769    
## as.factor(year)1997  0.31613859  0.25408695  1.2442      0.215485    
## as.factor(year)1998  0.46706610  0.27483535  1.6994      0.091441 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    33.616
## Residual Sum of Squares: 13.35
## R-Squared:      0.60286
## Adj. R-Squared: 0.49583
## F-statistic: 8.91819 on 24 and 141 DF, p-value: < 2.22e-16
###Fourth Model: Maximum Likelihood Estimator###
# For Comparison: Linear Model in MLE
linear.mle <- function(par, X, y, n) {
    beta<-par[-length(par)]
    sigma2<-par[length(par)]
    loglikelihood <- -(n/2)*log(2*pi)-(n/2)*log(sigma2)-
      (t(y-X%*%beta)%*%(y-X%*%beta)/(2*sigma2))
    return(loglikelihood)
}


#subset data in matrix form
X.mod<-as.matrix(subset(fh,select=c(constant,lnlmtue_1,DENSITY,DEIND,lngdp_pc,
                      UR,TRADE,FDI,LLVOTE,LEFTC,TCDEMC,GOVCON,OLDAGE,
                      y88,y89,y90,y91,y92,y93,y94,y95,y96,y97,y98)))

# do a linear model with a time lag and time dummies
lin.mod <- optim(par=c(0,0,0,0,0,0,0,0,0,0,0,0,
                       0,0,0,0,0,0,0,0,0,0,0,0,1),   # starting value for beta and sigma
   linear.mle,                  # the log-likelihood function
   method="BFGS",                   # optimization method
   hessian=TRUE,                    # return numerical Hessian
   control=list(fnscale=-1),        # maximize function instead of minimize
   y=fh$lnlmtue, X=X.mod, n=180)    # the data
## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced
# standard error of estimator
OI<-solve(-1*lin.mod$hessian)
se<-sqrt(diag(OI))
lin.mod$par; se
##  [1]  1.53859740656  0.75711795809  0.00075904866  0.00415123634
##  [5] -0.01824582781 -0.02507694984 -0.00003919987 -0.01409629595
##  [9]  0.00412294624  0.00095393272  0.00096184546  0.02718519240
## [13] -0.03505748214  0.09555748237  0.27821408818  0.11797066782
## [17] -0.09462528441  0.00215996393 -0.06591424593  0.05600699917
## [21]  0.13280965776  0.08864408076  0.14559411510  0.26519338209
## [25]  0.08645490262
##  [1] 0.7790569614 0.0489823918 0.0019540951 0.0135840395 0.1859839301
##  [6] 0.0111310886 0.0016638735 0.0190482414 0.0098078452 0.0006761811
## [11] 0.0012633199 0.0112462704 0.0222351484 0.1079039004 0.1098943158
## [16] 0.1126536122 0.1126779789 0.1116047453 0.1138235388 0.1146197103
## [21] 0.1154613847 0.1186928176 0.1219921841 0.1395409277 0.0091100291
#compare OLS to MLE
lin.mod.ols<-lm(lnlmtue~lnlmtue_1+DENSITY+DEIND+lngdp_pc+UR+TRADE+FDI+
                  LLVOTE+LEFTC+TCDEMC+GOVCON+OLDAGE+as.factor(year),data=fh)
summary(lin.mod.ols)
## 
## Call:
## lm(formula = lnlmtue ~ lnlmtue_1 + DENSITY + DEIND + lngdp_pc + 
##     UR + TRADE + FDI + LLVOTE + LEFTC + TCDEMC + GOVCON + OLDAGE + 
##     as.factor(year), data = fh)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.27217 -0.14718  0.01391  0.14856  1.19934 
## 
## Coefficients:
##                        Estimate  Std. Error t value Pr(>|t|)    
## (Intercept)          1.53815566  0.83681661   1.838   0.0679 .  
## lnlmtue_1            0.75715953  0.05261397  14.391   <2e-16 ***
## DENSITY              0.00075865  0.00209897   0.361   0.7183    
## DEIND                0.00415689  0.01459117   0.285   0.7761    
## lngdp_pc            -0.01822845  0.19977287  -0.091   0.9274    
## UR                  -0.02508193  0.01195635  -2.098   0.0375 *  
## TRADE               -0.00003922  0.00178723  -0.022   0.9825    
## FDI                 -0.01409620  0.02046049  -0.689   0.4919    
## LLVOTE               0.00412575  0.01053500   0.392   0.6959    
## LEFTC                0.00095433  0.00072631   1.314   0.1908    
## TCDEMC               0.00096110  0.00135698   0.708   0.4798    
## GOVCON               0.02717781  0.01208007   2.250   0.0259 *  
## OLDAGE              -0.03506607  0.02388367  -1.468   0.1441    
## as.factor(year)1988  0.09556586  0.11590395   0.825   0.4109    
## as.factor(year)1989  0.27819170  0.11804193   2.357   0.0197 *  
## as.factor(year)1990  0.11791994  0.12100580   0.974   0.3313    
## as.factor(year)1991 -0.09464141  0.12103198  -0.782   0.4354    
## as.factor(year)1992  0.00218939  0.11987917   0.018   0.9855    
## as.factor(year)1993 -0.06588254  0.12226247  -0.539   0.5908    
## as.factor(year)1994  0.05602039  0.12311767   0.455   0.6497    
## as.factor(year)1995  0.13276116  0.12402175   1.070   0.2861    
## as.factor(year)1996  0.08858347  0.12749276   0.695   0.4882    
## as.factor(year)1997  0.14550139  0.13103674   1.110   0.2685    
## as.factor(year)1998  0.26512769  0.14988656   1.769   0.0789 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3158 on 156 degrees of freedom
## Multiple R-squared:  0.8863, Adjusted R-squared:  0.8696 
## F-statistic: 52.89 on 23 and 156 DF,  p-value: < 2.2e-16
## Log-Likelihood Fuction for Spatial and Time Lag##
spatiotemporal <- function(par, X, y, W, T, N) {
    n<-N*T
    rho<-par[1]
    sigma2<-par[2]
    beta<-par[3:length(par)]
    ident<-diag(1,nrow=T)
    nIdent<-diag(1,nrow=N)
    bigIdent<-ident%x%nIdent
    bigW<-(ident%x%W)
    first<-prod(1-rho*eigen(bigW)$values) #Similar values for bigW or W. NaNs produced, though.
#   first<-det(bigIdent-rho*bigW)
    loglikelihood <- log(first)-(n/2)*log(2*pi)-(n/2)*log(sigma2)-
      ((t(y-rho*bigW%*%y-X%*%beta)%*%(y-rho*bigW%*%y-X%*%beta))/(2*sigma2))
    return(loglikelihood)
}

#prep data for simple model
inds<-model.matrix(~factor(fh$cc)-1)
X.mod.b.0<-as.matrix(subset(fh,select=c(constant,lnlmtue_1,y88,y89,y90,
                                        y91,y92,y93,y94,y95,y96,y97,y98)))
X.mod.b<-cbind(X.mod.b.0,inds)

#estimate simple model
spatiotemp.mod <- optim(par=c(0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
                              0,0,0,0,0,0,0,0,0,0,0),   # starting value for beta and sigma
   spatiotemporal,                  # the log-likelihood function
   method="BFGS",                   # optimization method
   hessian=TRUE,                    # return numerical Hessian
   control=list(fnscale=-1),        # maximize function instead of minimize
   y=fh$lnlmtue, X=X.mod.b,W=wmat, N=15, T=12)   
## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced
spatiotemp.mod$par
##  [1] -0.25897386  0.07752326  3.88603226  0.65694251  0.13879425
##  [6]  0.39534813  0.28949429  0.03774227  0.08594637 -0.06105515
## [11]  0.03530374  0.11915186  0.07396601  0.13881230  0.26196905
## [16]  0.41322640  0.33317484  1.03345956  0.68676933  0.30464497
## [21]  0.57479921 -1.87122859 -0.03187805  0.55175184  0.74181189
## [26] -0.07687160 -0.11597257  0.90109413  0.35857262  0.08267825
# standard error of estimator
OI<-solve(-1*spatiotemp.mod$hessian)
se<-sqrt(diag(OI))
## Warning in sqrt(diag(OI)): NaNs produced
spatiotemp.mod$par; se
##  [1] -0.25897386  0.07752326  3.88603226  0.65694251  0.13879425
##  [6]  0.39534813  0.28949429  0.03774227  0.08594637 -0.06105515
## [11]  0.03530374  0.11915186  0.07396601  0.13881230  0.26196905
## [16]  0.41322640  0.33317484  1.03345956  0.68676933  0.30464497
## [21]  0.57479921 -1.87122859 -0.03187805  0.55175184  0.74181189
## [26] -0.07687160 -0.11597257  0.90109413  0.35857262  0.08267825
##  [1] 0.067119427 0.008221808         NaN 0.055020379 0.102583244
##  [6] 0.107165680 0.113019875 0.110350365 0.106780185 0.104265401
## [11] 0.102911500 0.102835229 0.103713089 0.104187654 0.105700740
## [16]         NaN         NaN         NaN         NaN         NaN
## [21]         NaN         NaN         NaN         NaN         NaN
## [26]         NaN         NaN         NaN         NaN         NaN
##prep data for fuller model
inds<-model.matrix(~factor(fh$cc)-1)
X.mod<-as.matrix(subset(fh,select=c(constant,lnlmtue_1,DENSITY,DEIND,
                      lngdp_pc,UR,TRADE,FDI,LLVOTE,LEFTC,TCDEMC,GOVCON,
                      OLDAGE,y88,y89,y90,y91,y92,y93,y94,y95,y96,y97,y98)))
X.mod.a<-cbind(X.mod,inds)

#fit simple model
spatiotemp.full <- optim(par=c(0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
                               0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
                               0,0,0,0,0,0,0,0,0),   # starting value for beta and sigma
   spatiotemporal,                  # the log-likelihood function
   method="BFGS",                   # optimization method
   hessian=TRUE,                    # return numerical Hessian
   control=list(fnscale=-1),        # maximize function instead of minimize
   y=fh$lnlmtue, X=X.mod.a,W=wmat, N=15, T=12)  
## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced

## Warning in log(sigma2): NaNs produced
spatiotemp.full$par
##  [1] -0.28958562211  0.06910988208  3.90211984452  0.48481701150
##  [5]  0.01117189971  0.04760042241 -0.97604626466 -0.09466563444
##  [9]  0.00025650397 -0.02996796219 -0.00155310619  0.00002044894
## [13] -0.00275479490  0.02704489833  0.03058603721  0.18134892304
## [17]  0.46678520401  0.39136404912  0.15000329490  0.21378639844
## [21]  0.14111368237  0.23186297483  0.29601755581  0.29336485560
## [25]  0.36276560985  0.53421634167  0.41002618666  0.30536282955
## [29]  0.75406967478  0.57026132394  0.70299558387  1.11431934970
## [33] -2.14164687950  0.22831008994  0.53968654931  0.52913420631
## [37] -0.53846008443  0.44463817520  0.30271805641  0.77162464170
## [41] -0.09091985909
# standard error of estimator
# insoluable--bootstrapping may be the best solution here

0.0.9 Spatial Point Pattern Analysis

#R code to test complete spatial randomness
#Source: Chapter 7 of Bivand, Pebesma, & Gomez-Rubio 2008
#Usage of this code requires citation of the original source

#clean up and set directory
rm(list=ls())

#load libraries
library(spatstat)
library(maptools)
library(rgdal)
library(maps)

#load three data sets
data(cells)
data(japanesepines)
data(redwood)

###Fun with Plotting###
#redefine as SpatialPoints objects
spcells<-as(cells,"SpatialPoints"); summary(spcells)
## Object of class SpatialPoints
## Coordinates:
##      min max
## [1,]   0   1
## [2,]   0   1
## Is projected: NA 
## proj4string : [NA]
## Number of points: 42
spjpines<-as(japanesepines,"SpatialPoints"); summary(spjpines)
## Object of class SpatialPoints
## Coordinates:
##      min max
## [1,]   0 5.7
## [2,]   0 5.7
## Is projected: NA 
## proj4string : [NA]
## Number of points: 65
spredwood<-as(redwood,"SpatialPoints"); summary(spredwood)
## Object of class SpatialPoints
## Coordinates:
##      min max
## [1,]   0   1
## [2,]  -1   0
## Is projected: NA 
## proj4string : [NA]
## Number of points: 62
#convert to a unit square
spcells1<-elide(spcells,scale=TRUE,unitsq=TRUE); summary(spcells1)
## Object of class SpatialPoints
## Coordinates:
##      min max
## [1,]   0   1
## [2,]   0   1
## Is projected: NA 
## proj4string : [NA]
## Number of points: 42
spjpines1<-elide(spjpines,scale=TRUE,unitsq=TRUE); summary(spjpines1)
## Object of class SpatialPoints
## Coordinates:
##      min max
## [1,]   0   1
## [2,]   0   1
## Is projected: NA 
## proj4string : [NA]
## Number of points: 65
spredwood1<-elide(spredwood,shift=c(0,1),scale=TRUE,unitsq=TRUE); summary(spredwood1)
## Object of class SpatialPoints
## Coordinates:
##      min max
## [1,]   0   1
## [2,]   0   1
## Is projected: NA 
## proj4string : [NA]
## Number of points: 62
#plot the three data sets from Figure 7.1 of Bivand, Pebesma, and Gomez-Rubio 2008
pdf('cells.pdf');plot(spcells1); axis(1); axis(2); box();dev.off()
## quartz_off_screen 
##                 2
pdf('pines.pdf');plot(spjpines1); axis(1); axis(2); box();dev.off()
## quartz_off_screen 
##                 2
pdf('redwood.pdf');plot(spredwood1); axis(1); axis(2); box();dev.off()
## quartz_off_screen 
##                 2
###testing complete spatial randomness (Note:the code in the "plotting" section is unnecessary for this)###
#For the Japanese Pines
envjap<-spatstat::envelope(japanesepines,fun=Gest,nrank=2,nsim=99) #gets confused relative to "boot" library's version
## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.
#plot(envjap)
plot(y=envjap$obs,x=envjap$theo,type='l')
lines(y=envjap$lo,x=envjap$theo,lty=3,col='red')
lines(y=envjap$hi,x=envjap$theo,lty=3,col='red')

##Use the F function: distance from a point to the nearest event##
Fenvjap<-spatstat::envelope(japanesepines,fun=Fest,nrank=2,nsim=99)
## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.
plot(y=Fenvjap$obs,x=Fenvjap$theo,type='l')
lines(y=Fenvjap$lo,x=Fenvjap$theo,lty=3,col='red')
lines(y=Fenvjap$hi,x=Fenvjap$theo,lty=3,col='red')

###Redo for Redwoods###
envred<-spatstat::envelope(redwood,fun=Gest,nrank=2,nsim=99)
## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.
plot(y=envred$obs,x=envred$theo,type='l')
lines(y=envred$lo,x=envred$theo,lty=3,col='red')
lines(y=envred$hi,x=envred$theo,lty=3,col='red')

Fenvred<-spatstat::envelope(redwood,fun=Fest,nrank=2,nsim=99)
## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.
plot(y=Fenvred$obs,x=Fenvred$theo,type='l')
lines(y=Fenvred$lo,x=Fenvred$theo,lty=3,col='red')
lines(y=Fenvred$hi,x=Fenvred$theo,lty=3,col='red')

###Redo for Cells###
envcell<-spatstat::envelope(cells,fun=Gest,nrank=2,nsim=99)
## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.
plot(y=envcell$obs,x=envcell$theo,type='l')
lines(y=envcell$lo,x=envcell$theo,lty=3,col='red')
lines(y=envcell$hi,x=envcell$theo,lty=3,col='red')

Fenvcell<-spatstat::envelope(cells,fun=Fest,nrank=2,nsim=99)
## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.
plot(y=Fenvcell$obs,x=Fenvcell$theo,type='l')
lines(y=Fenvcell$lo,x=Fenvcell$theo,lty=3,col='red')
lines(y=Fenvcell$hi,x=Fenvcell$theo,lty=3,col='red')

###############################################################

#Source: Monogan, Konisky, & Woods (N.d.)
#Usage of the code or data requires citation of the original source


############## ESTIMATING MODELS #########################
#clean up and set directory
rm(list=ls())
options(scipen=10,digits=10)

#load libraries
library(spatstat)
library(maptools)
library(rgdal)
library(maps)
library(foreign)
library(mgcv)
library(rgeos)
#library(xtable)
library(akima) #for drawing image/surface plots
library(scatterplot3d)
#library(spdep)
checkCRSArgs <- function(uprojargs) {
  .Call("checkCRSArgs", uprojargs, PACKAGE="rgdal")
  }

#Load data
full<-read.csv('mnAir.csv')

#GAM model
gam.air<-gam(air~distances+s(longitude,latitude), data=full, family=binomial(link='logit')); summary(gam.air)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## air ~ distances + s(longitude, latitude)
## 
## Parametric coefficients:
##                  Estimate    Std. Error  z value   Pr(>|z|)    
## (Intercept)  1.0941008720  0.3020727202  3.62198 0.00029236 ***
## distances   -0.0021650171  0.0009046893 -2.39311 0.01670643 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                          edf   Ref.df  Chi.sq p-value
## s(longitude,latitude) 4.3121 5.969084 9.31141 0.15312
## 
## R-sq.(adj) =  0.0318   Deviance explained = 3.11%
## UBRE = 0.32832  Scale est. = 1         n = 548
for.comparison<-gam(air~distances+s(eastings,northings), data=full, family=binomial(link='logit')); summary(for.comparison)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## air ~ distances + s(eastings, northings)
## 
## Parametric coefficients:
##                 Estimate   Std. Error  z value   Pr(>|z|)    
## (Intercept)  1.100839924  0.301218429  3.65462 0.00025756 ***
## distances   -0.002185183  0.000901764 -2.42323 0.01538311 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                            edf   Ref.df  Chi.sq p-value
## s(eastings,northings) 4.435246 6.112045 9.52867 0.15314
## 
## R-sq.(adj) =  0.0325   Deviance explained = 3.18%
## UBRE = 0.32787  Scale est. = 1         n = 548
#draw figure of smoothed intercept
mn<-map("state","minnesota",resolution=0)

vis.gam(gam.air,view=c("longitude","latitude"),color='gray',
        plot.type='contour',type='link',cond=list(distance=0),main="", 
        xlab="Longitude", ylab="Latitude",ylim=mn$range[3:4],xlim=mn$range[1:2])
lines(mn,col='red',lwd=2)
points(x=full$longitude[full$air==1], y=full$latitude[full$air==1], pch='+',col='white')
points(x=full$longitude[full$air==0], y=full$latitude[full$air==0], pch=1,col='white')

############# Testing for complete spatial randomness. If it doesn't fail, something is bizarre. ###################

#Draw a map of the 50 states
state<-data(stateMapEnv)
us.states<-map("state", fill=TRUE, plot=TRUE, resolution=0) 

IDs <- sapply(strsplit(us.states$names, ":"), function(x) x[1])
us.map.0<-map2SpatialPolygons(us.states, IDs=IDs,
                              proj4string=CRS("+proj=longlat +datum=NAD27"))
us.map<-spTransform(us.map.0, CRS("+proj=longlat +datum=NAD83")) 

#extract minnesota
mn<-SpatialPolygons(Srl=list(us.map@polygons[[22]]),
                    proj4string=CRS("+proj=longlat +datum=NAD83"))
reach<-mn@bbox[1,2]

#prepare for CSR test
mn.trim<-mn@polygons[[1]]@Polygons[[1]]@coords
mn.coords<-dim(mn.trim)[1]-1
mn.bound<-owin(xrange=c(min(mn.trim[,1]),max(mn.trim[,1])),
               yrange=c(min(mn.trim[,2]),max(mn.trim[,2])),
               poly=list(x=rev(mn.trim[1:mn.coords,1]),
                         y=rev(mn.trim[1:mn.coords,2])))
mn.ppp<-ppp(x=full$longitude, y=full$latitude, window=mn.bound)
## Warning: 9 points were rejected as lying outside the specified window
## Warning: data contain duplicated points
# There are 548 observations, but 8 are tossed for being illegal.

#look at it
plot(mn.ppp)
## Warning in plot.ppp(mn.ppp): 9 illegal points also plotted

#conduct CSR tests
env.mn<-spatstat::envelope(mn.ppp,fun=Gest,nrank=2,nsim=99); env.mn
## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.
## Pointwise critical envelopes for G(r)
## and observed value for 'mn.ppp'
## Edge correction: "km"
## Obtained from 99 simulations of CSR
## Alternative: two.sided
## Significance level of pointwise Monte Carlo test: 4/100 = 0.04
## .....................................................................
##      Math.label     Description                                      
## r    r              distance argument r                              
## obs  hat(G)[obs](r) observed value of G(r) for data pattern          
## theo G[theo](r)     theoretical value of G(r) for CSR                
## lo   hat(G)[lo](r)  lower pointwise envelope of G(r) from simulations
## hi   hat(G)[hi](r)  upper pointwise envelope of G(r) from simulations
## .....................................................................
## Default plot formula:  .~r
## where "." stands for 'obs', 'theo', 'hi', 'lo'
## Columns 'lo' and 'hi' will be plotted as shading (by default)
## Recommended range of argument r: [0, 0.15055]
## Available range of argument r: [0, 0.41665]
#env.mn<-spatstat::envelope(mn.ppp,fun=Gest,nrank=2,nsim=99,savepatterns=TRUE); env.mn
plot(y=env.mn$obs,x=env.mn$theo,type='l',xlim=c(0,1),ylim=c(0,1))
lines(y=env.mn$lo,x=env.mn$theo,lty=3,col='red')
lines(y=env.mn$hi,x=env.mn$theo,lty=3,col='red')

Fenv.mn<-spatstat::envelope(mn.ppp,fun=Fest,nrank=2,nsim=99); Fenv.mn
## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.
## Pointwise critical envelopes for F(r)
## and observed value for 'mn.ppp'
## Edge correction: "km"
## Obtained from 99 simulations of CSR
## Alternative: two.sided
## Significance level of pointwise Monte Carlo test: 4/100 = 0.04
## .....................................................................
##      Math.label     Description                                      
## r    r              distance argument r                              
## obs  hat(F)[obs](r) observed value of F(r) for data pattern          
## theo F[theo](r)     theoretical value of F(r) for CSR                
## lo   hat(F)[lo](r)  lower pointwise envelope of F(r) from simulations
## hi   hat(F)[hi](r)  upper pointwise envelope of F(r) from simulations
## .....................................................................
## Default plot formula:  .~r
## where "." stands for 'obs', 'theo', 'hi', 'lo'
## Columns 'lo' and 'hi' will be plotted as shading (by default)
## Recommended range of argument r: [0, 0.42564]
## Available range of argument r: [0, 0.42564]
plot(y=Fenv.mn$obs,x=Fenv.mn$theo,type='l',xlim=c(0,1),ylim=c(0,1))
lines(y=Fenv.mn$lo,x=Fenv.mn$theo,lty=3,col='red')
lines(y=Fenv.mn$hi,x=Fenv.mn$theo,lty=3,col='red')

#Very clear patterns of clustering on both counts.

#Fit a point process model
fit.mn<-ppm(mn.ppp,~1); summary(fit.mn)
## Point process model
## Fitting method: maximum likelihood
## Model was fitted analytically
## Call:
## ppm.ppp(Q = mn.ppp, trend = ~1)
## Edge correction: "border"
##  [border correction distance r = 0 ]
## ---------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
## 
## Data pattern:
## Planar point pattern:  539 points
## Average intensity 21.1 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 353 vertices
## enclosing rectangle: [-97.2252121, -89.47309113] x [43.49322891, 
## 49.38323212] units
## Window area = 25.5331 square units
## Fraction of frame area: 0.559
## 
## Dummy quadrature points:
##      64 x 64 grid of dummy points, plus 4 corner points
##      dummy spacing: 0.12112689018 x 0.09203130007 units
## 
## Original dummy parameters: =
## Planar point pattern:  2367 points
## Average intensity 92.7 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 353 vertices
## enclosing rectangle: [-97.2252121, -89.47309113] x [43.49322891, 
## 49.38323212] units
## Window area = 25.5331 square units
## Fraction of frame area: 0.559
## Quadrature weights:
##      (counting weights based on 64 x 64 array of rectangular tiles)
## All weights:
##  range: [0.000429, 0.0111]   total: 25.5
## Weights on data points:
##  range: [0.000429, 0.00557]  total: 1.47
## Weights on dummy points:
##  range: [0.000429, 0.0111]   total: 24
## ---------------------------------------------------------------------------
## FITTED MODEL:
## 
## Stationary Poisson process
## 
## ---- Intensity: ----
## 
## 
## Uniform intensity:
## [1] 21.10985456
## 
##                Estimate          S.E.     CI95.lo     CI95.hi Ztest
## log(lambda) 3.049739972 0.04307304923 2.965318347 3.134161597   ***
##                    Zval
## log(lambda) 70.80390237
## 
## ----------- gory details -----
## 
## Fitted regular parameters (theta):
## log(lambda) 
## 3.049739972 
## 
## Fitted exp(theta):
## log(lambda) 
## 21.10985456
plot(simulate(fit.mn,2))
## Generating 2 simulated patterns ...1,  2.

#beautiful, the simulation (also used in "envelope") respects the polygon boundary